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Two methods of simulation of multi correlated random processes 
from the given matrix of spectral density function have been pre- 
sented. It has been noted that FFT method works as efficiently as 
the Trigonometric method and is much faster. It has been found that 
there are certain cases in which Trigonometric approach has advan- 
tages over the FFT method. Some example problems are solved to show 
the usefulness of this approach in solving the problems of linear 
and nonlinear random vibrations. It has been observed that this 
technique and particularly the FFT method offers a very fast and 
convenient alternative for performing random nonlinear response 
analysis. Various possible areas in which this approach can be 
extended have been also discussed. 
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CHAPTER I 


INTRODUCTION 

In this dissertation, two methods are described to simulate on a 
digital computer a set of correlated, stationary and Gaussian time 
series with zero mean from the given matrix of power spectral densi- 
ties and cross-spectral densities. Some example problems are investi- 
gated to show the power of this technique to solve the problems of 
linear and in particular nonlinear random vibrations. The development 
set forth is for any arbitrary number of correlated series; however, 
in practice, the number is limited by the storage capability of the 
computer. 

The first method is based upon trigonometric series with random 
amplitudes and deterministic phase angles. The random amplitudes are 
generated by using a standard random number generator subroutine. An 
example is presented which corresponds to three components of wind 
velocities at two different spatial locations for a total of six 
correlated time series. Selected spectral densities computed from the 
simulated time series are compared to the original spectral densities 
from which the time series were generated. 

In the second method, the whole process is carried out using the 
Fast Fourrier Transform approach in place of trigonometric series. It 
is found that this method gives more accurate results and works about 
twenty times faster for a set of six correlated time series. 
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To show one of the many areas of application of the present 
method of simulation, namely the class of random structural vibration 
analysis, the following problems have been investigated: (1) The 

linear vibration characteristics of a long tower under the action of 
correlated wind loads have been presented. Taking it to be a fourteen 
degrees of freedom system, the time history of the displacement of the 
top of the tower has been plotted considering up to three modes of 
vibration. The usefulness of this method in the case of linear 
analysis can be significant and very important, e.g., looking for the 
occurrence of maximum structural response rather than the r.m.s. value 
of response. This knowledge will be very useful for the reliability 
study of structures under random loads; (2) One of the most interesting 
and significant applications of the proposed method is the simulation 
of random generalized forces. The necessity of simulating random 
generalized forces arises when the dynamic response analysis is per- 
formed in time domain either for the purpose of obtaining information 
beyond the second order statistics (such as first passage time distri- 
bution) or when the structure is nonlinear and therefore an approximate 
random response is sought by simulating the excitation. 

In order to assert the validity of the preceding discussion, the 
problem of the nonlinear vibration of a string has been solved in the 
time domain by simulating the random generalized forces. 

Finally to show the application of this method to a more complex 
problem, the random nonlinear vibration of a flexible plate immersed 
in a fluid flow on one side and backed by a fluid filled cavity of 



3 


finite dimensions on the other side is considered. The nonlinear 
plate stiffness induced in the plate by out-of-plane bending and the 
mutual interaction between the external and the internal fluid flow 
is included. 

The FFT simulation technique is utilized for the response 
analysis of the plate undergoing large deformation. The same problem 
has been done by Shinozuka [1] where he has taken a multidimensional 
trigonometric series model for the generalized forces. The analysis 
is performed in the time domain rather than in the frequency or wave 
number domain as is usually done in linear response analysis. The 
numerical example has been presented for subsonic flow over the plate. 



CHAPTER II 


LITERATURE SURVEY 

Part A : Review of Existing Simulation t’lethods 

Numerous papers dealing with simulation of random process have 
been published in recent years. Although many authors dealt with the 
simulation of single random processes utilizing trigonometric series 
[2], filtered white noise [2], filtered shot noise [3] and correlated 
random pulse trains [4] etc, only Hoshiya and Tideman [5], Shinozuka [6] 
and Borgman [7] studied the simulation of multi correlated processes. 

I. In the simulation of ocean surface elevation, Borgman used wave 
superposition by choosing the frequency in such a way that the ampli- 
tude of each wave function was an equal portion of the cumulative 
spectrum. Borgman also presented a method for simulating several 
simultaneous time series by passing a white noise vector through 
filters. He proposed the following model [7]. 

m 

«»00 

m = 1, 2, ... M (2-1) 


where. 
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x^(t) = independent random inputs 
kjjjj(-t) = kernals 

Let represent the cross spectral density between y^(t) and 

yj-(t). Kernals kj^j (t) and its Fourier Transform Kj^j (f ) are deter- 
mined from the relation 

S^r(f) = r Kj„j(f)K^j(f)S^.(f), 

3=1 

r = 1 , 2, , . . , m and (2-2) 

m = 1 , 2, . . . , M 

where the bar denotes the complex conjugate. 

This system of equation can be solved sequentially (taking Kjj to 
have zero gain for j = 1 , 2, 3 . . . ) as 
Kn(f) = v^TTF) 

K2i(f) = /S2l(^') (2-3) 

K22(f) = [522(f) - lK2i(f)|^]^*etc. 

This, in turn, determines the digital filter co-efficients, a^p^j, 
needed to approximate the kernel associated with the system function 
kmj(f). Then the final simulation equation reduces to 

N N 

y„(kAt) = E I a ^Xj, k-n, 
j=l n=-N 

m = 1 , 2, ... M (2-4) 

in which Xj, n for n = 1 , 2, 3, ... is the j th generated sequence of 
independent, zero mean, unit variance, normal random variables. 
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II. Later Hoshiya and Tieleman [5] considered the following trigo- 
nometric model for simulating two correlated random processes: 

n 


x(t) = 


E 

j=l 


(ajcos Wjt + bjSin Wjt) 


n 

y(t) = E 
3=1 


{CjCOS (Wjt + aj) + djSin (Wjt + aj) 


(2-5) 


where aj, bj , Cj and dj are random variables and aj is deterministic. 
The two stationary normal processes are correlated if there exists a 
non-zero correlation between aj and Cj and/or bj and dj. The co- 
variance Cj^yj between the two random variables aj and Cj is defined 


^xyj "" 

and the standard deviation at frequency Wj, a^j and Oyj are given as 
Oxj = ’^2S^(wj )Aw (2-6) 

0yj = ’^2Sj(Wj)Aw (2-7) 

where S^(Wj) and Sy(wj) are the discrete form of the power spectra for 
processes x and y respectively. 


Expressions for the co-spectrum and the quadrature spectrum are given as 


2Q^y(wj)Aw = ■<^xj‘^yjPxyj5'*'^ “j 


whence. 


-1 

aj = tan 


Cxy(Wj)] 


(2-8) 

(2-9) 


( 2 - 10 ) 


and 

Pxyj = {C^y(w.j)} {Qxy(w)}^ 

S$(wj)Sj(wj) 


( 2 - 11 ) 
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The random variables a. and are generated independently from a 
normal distribution with zero mean and a standard deviation 

= ’^2Sx(w)Aw and = *^2S^(wj)Aw respectively. Then Cj and dj 
are generated as follows. Since aj and Cj are both normally dis- 


tributed, the joint probability density function of a^ and c^ is 

/ . r ...2 ^ 

P(aj.Cj) = exp 


2TTaxjOyj 


^Pxyj^j^.i ^ fj 

^xj^yj ^yj^ 


{* 

]} 


12. 


^xj 


(2-lla) 


and the probability density function of a,- is 

r =2 - 


I UIIUU lUli U I 

{-*} 


(2-llb) 


The conditional probability density function of Cj for given aj is 


P(c./aj) = = 

■ '-J 




X exp 


(■ 


2a 


yj 


’^^■Pxyj 




XJ 


(2-llc) 


Consequently, the conditional probability density function of c. for 

Oyi 

given aj is also Gaussian with mean Pxyj^'^j ® standard deviation 

XJ 


c^u-i’^l-Pyvi • Similarly, the conditional distribution of d- for a given 
yj ^yj o • 

b- is Gaussian with mean Pxyj^^j and standard deviation of c^j*^l -p^y?. 

XJ 

Thus random variable Cj is generated from Gaussian distribution with 


mean 


Pxyj^^j and a standard deviation of Pyj*^-Pxyj* Similarly, the 

X J 

random variable dj is generated from a Gaussian distribution with mean 
Pxyj^^j and a standard deviation of Pyj*^"Pxyj' 
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III. In his paper [8], Shinozuka proposed a different 
trigonometric model for the simulation of multivariate random pro- 
cesses. He used a series of cosine functions with weighted ampli- 
tudes, almost evenly spaced random frequencies and random phase 
angles with uniform distribution. He considered a set of n sta- 
tionary random processes fo^(t) (i = 1, 2, ... n) with a specified 
cross-spectral density matrix S°(w) = [S°f.f.(w)l, where S°f.^.(w) 

are mean square spectral densities if i = 3 and cross spectral 
densities of fo.j(t) and foj(t) if i 4= j. 

He proposed the following model 


i , 2^15 N 


X cos 

(2-12) 

where, 


Wji^ (k = 1, 2, ... N) are random variables identically and independent 

ly distributed with the density function 


gj(w) = |Hjj(w)l^/aj^ 

00 

/ 

(2-13) 

with = J lHjj(w)|^dw 

-00 

(2-14) 

and Qji^ (k = 1, 2, ... N) are also identically and 

independently dis- 


tributed with the uniform density 1/(27 t) between 0 and 2 tt . Hj^(w) 
are obtained from 

SijM = 5 

K“ I 


(2-15) 
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and. 


Yij(w) = 

e,j(w) = t 


lReH,7(i?)l ■ 


(2-16) 

(2-17) 


The series simulated by the above technique have the asymptotic 
ergodicity. 

IV. Comparison of Various Approaches . 

It appears that method of simulation proposed by Shinozuka [8] is 
the most efficient so far. Thismodel [8] seems to be more general than 
that due to Hoshiya and is more efficient than that due to Borgman [7]. 
Since the measured cross spectral densities are usually given 
numerically in terms of the real and the imaginary parts or the 
moduli i and the angles, the computation work associated with finding 
lHij(w)l, from which 7 -jj(w) = iH-j j(w) |/Hjj (w) and the angles 9ij(w) 
can be directly evaluated in the case of Ref. [8]. Therefore, method 
of simulation presented here appears more practical than that proposed 
in Ref. [7], which requires (a) the inverse Fourier Transformation of 
N(N+1)/2 functions of w, Hij(w) and (b) the same number of integra- 
tion in the time domain. Also, form of simulated functions con- 
sisting of sum of cosine functions can be proved extremely advan- 
tageous, when used as input, in evaluating the corresponding output 
of a linear system. 
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PART B . 

The main usefulness of the simulation method is in the area of 

(a) a numerical analysis of dynamic response of non-linear structures, 

(b) time domain analysis of linear structures under random excitations 
performed for the purpose of obtaining the kind of information, such as 
first excursion probability and exact time history of a sample 
function, that is not obtainable from the standard frequency domain 
analysis, and (c) numerical solution of stress wave propagation 
through a random medium and eigen value problems of structures with 
randomly non-homogeneous material properties. The usefulness of this 
method has been demonstrated in the later chapter of this dissertation. 
In this connection, it looks appropriate to make a relative study of 
the prevalent methods of handling the non-linear random vibration 
problems. 

Review of Stationary Random Response of 
Multi degree of Freedom Nonlinear Systems 

Fokker-Planck Approach . 

One exact method of studying the stationary random response of 
a nonlinear system is the Fokker-Planck approach. If the excitation 
is a Gaussian white noise, then the transitional probability density 
Of the response process is governed by the Fokker-Planck equation. 

The equation governing the first probability density function for the 
stationary response has been solved under the following rather 
restrictive conditions: 
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(1) the damping force is proportional to the velocity 

(2) the excitation is a Gaussian white noise 

(3) the correlation function matrix of the excitation is 

proportional to the damping matrix of the system. Under the above 

conditions, the equation of motion may be written as follows: 

Mx + Cx + = f(t) (2-18) 

3x 

wi th 



ny = 0 

R:p<T) = 2 yC6(t) (2-19a) 

where y is a constant, u(7) is the potential energy of the system and 

8u 

(2-19b) 
3x 

Suppose that there exists an orthogonal matrix A which can simul- 
taneously diagonal ize the mass matrix M and the damping matrix C; 
a''’a = I ^ 

A% = V 1 (2-20) 

A^CA = A ) 

where V and A are two diagonal matrices. Then, upon using the trans- 
formation X = Az and noting that 


3x 8x 3z 


Mil 

m = 1 , ... n 
^^m 

= 3ujz) . 

3z 


(2-21) 
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Equation (2-18) becomes 


3u(¥) t 

Vz + AZ + A'f(t) = b(t) 


(2-22a) 


and the correlation function matrix of b(t) is given by 

= 2yA<S(t) . (2-22b) 
The stationary Fokker Planck equation associated with (2-22a) is given 
by 


E 

j= 


a , , % a /r • 3u(i)]V 


-0 

- E 2 0 = 0 


(2-23) 


where X- and v- denote the j th diagonal element of A and V, p is the 
J j j 

abbreviation for the first probability density of the Markovian 


vector! 


{!)• 


The solution to (2-23) may be written as follows: 


p(z,z) = S exp< - - 




(2-24) 


This solution was first obtained by Ariarathan [9] for a two-degree- 
of freedom system, and it was extended to the above form by Caughey 
[10]. The constant 8 in (2-24) is a normalizing factor such that 


p(z, 7 )dzi . . .dZj,dz^ . . .dZj^ = !• 

In the original co-ordinates, equation (2-24) becomes 


(2-25) 


p(x,x) = B exp 


/- + u 

^ .1. J.I j • 


(30 


(2-26) 


It will be noted that the terms in the square brackets are 
respectively the kinetic energy and the potential energy of the system. 



13 


Equation (2-26) may also be written as 

p(x,^ = 6 exp xM^x} exp {- ^(x)} (2-27) 

Hence x and x are linearly independent. 

Normal Mode Approach 

Consider an n-degree-of-freedom system governed by the equation 
of motion 

+ K^o))T + y^()T,7) = f(t), (2-28) 

u = a small parameter 

The matrices and are respectively the damping matrix and 
the stiffness matrix of the system due to the linear part of the 
damping forces and the spring forces, and represents the non- 

linear forces of the system. Here, f(t) is a stationary Gaussian random 
vector. Without loss of generality, one can assume that the mean vector 
of f m^ = 0. 

In using this approach, the following two conditions must be 

satisfied: (1) the linear system obtained by neglecting the nonlinear 

* 

term ^()T,jr) in (2-28) must possess normal modes; (2) the correlation, 
function matrix R:p{T) must be diagonalized by the same matrix which 
diagonalizes the matrices M, and 

The second condition is quite restrictive and is seldom realized in 
real systems . 

Assume that the above restriction can be met, then there exists 


a matrix A such that 
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a'^’ma = I 
aTk^°U = 

a'*’c^o)a = A^o), Xkj^°^ = Xk(°)6y 
A^Rf(T)A = D(t), d|^j = d|^(T)6|(j . 

By using the transformation x = A7, Eq. (2-28) reduces to 
z" + A^°^7 + + yA^g(z,7) = A^T(t) = b(t) 

Where the correlation function matrix of F is 
R^t) = D(t) . 

In component form Eq. (2-30a) becomes 


(2-29) 


(2-30a) 

(2-30b) 


and Eq. (2-30b) becomes 

E[b|^(t)bj (t+x) ] = d|^(T)6j^j, j»k“l» ••• ri» 
The differential equation in (2-31 a) may be written as 

^ '•iS * 


where the deficiency term ej is given by 


j = 1. 

n _ 

+ y^E^a|^jg|^(z,z), j = 1 , . . . n 


(2-31a) 


(2-31b) 

(2-32) 


(2-33) 


2 

If the quantities Xj and Wj are chosen in such a way that some 
measure of the deficiency term is minimized, then it seems reasonable 
that the statistics of the response of the nonlinear system can be 
approximated by those of the linear system described by 


(2-34) 
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At this stage, the differential equations are uncoupled and the 
excitation b(t) is an uncorrelated vector process. Hence each un- 
coupled differential equation can be solved separately, 

o 

In order to determine Xj and Wj , Caughey [11] chose them so as 
to minimize the mean square value of the deficiency term e. This can 
be achieved by requiring that 



(2-35) 


Substituting (2-33) in (2-35) and interchanging the order of dif- 
ferentiation and expectation, we obtain 

+ y^Z^akjE[Zj-g,,(I.I)]/E[Zj2] 


= (w>))2 + u E ak^[ZjE[Zjg^(t,z)l (2-36) 

IC“ 1 

/E(z/] ■ 

Equations (2-34) and (2-36) can be used to find various mean square 
values of the response process. 

In certain cases, the contribution from the first mode may be 
dominant. In these cases, we may let Xj = aj-|Zi in the above 
derivation. Then 



(2-37) 
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This is rather a rough approximation, but it is very simple and in some 
cases, it does give reasonably approximate solutions. 

Perturbation Approach 

Consider the same problem as defined in the previous section, whose 
equations of motion are 

Mx + c(o)x + K^°^x + yg(x,x) = f(t) (2-38) 

Assume that y is so small that the solution of Eq. (2-38) can be 

approximately represented by [12] 

X = Xq + y)C| (2-39) 

2 

Substituting (2-39) into (2-38), neglecting terms involving y , 
y-^, ... and equating corresponding coefficients of y and y yields 
the following sets of linear differential equations: 

Mxq + c(°^Xg + k(°^Xq = f(t) (2-40a) 

Mx^ + C^°)x^ + k(o)x-| = -?(xo(t),Xo(t)) (2-40b) 

Correct to the same order of accuracy, the instantaneous correlation 
matrix for displacement becomes 

E[xx^] = E[xqXo^] + y{E[XQx/] + E[x^Xq^]}- (2-41) 

Note that 

E[XqX-|!*’]= (E[XqX-,^])^. (2-42) 

The matrix E[XqXq^] can be found from (2-40a) by the various approaches 
as in linear analysis, and E['XqX.|^] may be evaluated as follows. 

Since (2-40a) and C2-40b) are linear, their stationary solutions are 

Xq = G(t-T)f(i)dT 

-CO 


(2-43a) 
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fOO 


J 6(t-T)g(T)dT , 


(2-43b) 


where G(t) is the common impulse response function matrix of (2-40a) 
and (2-40b) and g(x) is the abbreviation of 9 (xq(t) ,Xq(t)) . Thus, 

/OO fO 


E[XqX-|^] = - 


G(t-T-|)E[f(T])g^(T2)] 


X G^(t-T2)dT-|dT2 • (2-44) 

The matrix E[f(T-| )^(t 2 )] can be evaluated with the help of proper- 
ties of Gaussian processes. Therefore, we can find E[^QXy^] and 
hence E[xx^]. 

Usually, the evaluation of the integrals in (2-44) is not easy, so 
Tung [13] has developed a different approach to generate E[xx*^] 
from (2-40a) and (2-40b). He applies Foss's method [14] to uncouple 
(2-40a) and (2-40b) into first order differential equations and then 
solves the resulting equations to find the various instantaneous 
correlation matrices. 

This approach will fail if the damping matrix is a null 
matrix. In this case. Equation (2-40a) does not have a stationary 
solution since all its correlation functions will finally go to 
infinity. Another limitation of this approach is that not only must the 
nonlinearity of the system has to be small, but also the excitation 
has to be sufficiently low. 

Generalized Equivalent Linearization 

The normal mode approach is quite powerful if it applies. How- 
ever, due to the conditions imposed on the excitation, its application 
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is rather limited. 

In his thesis Yang [15] introduced a more general approach. 
Except that the excitation must be stationary, the only additional 
restriction to this approach is thatthe excitation be Gaussian. In 
this approach, we define an auxiliary set of linear differential 
equations for the original nonlinear system. The solution of the 
original nonlinear system is approximated by the solution of an 
auxiliary set and the unknown co-efficients are chosen in such a way 
that some measure of the difference between the two sets of equations 
is a minimum. 

Consider the following linear equation 

Mx + Cx + Kx = f(t) (2-45) 


and the nonlinear system given by 
Mx + = T(t) . 

The difference between (2-46) and (2-45) 

• • 

? = g()T,)0 - cjT - Kx 

The necessary conditions to minimize Ele"^?] are given by 

9E[e^e] = 2E[e ^ 3e ] = 2E[e.x.] = 0 
“jk ‘ 5'^jk ' " 

3E[e^e] = 2Efe^ e 1 = 2E[e,x,.l = 0- 

‘ 3^jk f ^ 


(2-46) 

(2-47) 

(2-48) 


Upon using Eq. (2-47), they become in the matrix form 
E[ex^] = E[g(x,^^^] - E[xx^] - KE[xx^] = 0 
E[ex^l = E [^(x ,jT)>^^] - CE[xx^] - KE[xx^] = 0 
It has been shown that the conditions in Eq. (2-49) do define a 
minimum for E[¥^e][15]. 
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In order to solve Eq. (2-49) for K and C, it is first necessary to 
express E[g(x,x)x‘] and E[g(x,x)x 1 in terms of E[xx'], E[xx‘] and 
E[xx^]. Let denote the displacement of k th mass relative to 
the r th mass and let the approximate force acting on the k th mass 
by the nonlinear element connecting the k th mass and the r th mass 
be denoted by S|^y.(y,^^,y|^^) . Then 

E[g,^(x,x),Xj.] = E E[S|^^(y|^^,y|^^)Xj] 

(2-50) 

Elg|j( 7 . 7 )xjl = E E[S|^^(y|^^,y )xjl , 
r 

r^k 


where the sum is taken over all nonlinear elements connected to the 

k th mass. Since 7 is a Gaussian vector, it follows that the 

• « 

quantities yj^^, yj^^, Xj, and Xj will be Gaussian distributed. Hence 


E[S|^r(ykr’'^kr^^j^ " ^^^kr^'^kr’^kr^-^kr^ 

^ ••>'kr>^kr> ' ^ tJ'kr* j ' /^Cykrl 

^'\r'^kr->'kr>'‘j' ° V'^kr-J^kri^kr' 

X E[y|^^Xj]/E[y^^J 

+ E[S|^^(;k^.y|^^)y^^]-E[y^^ij]/E[y^^J 

Define 

• • * 2 

Tkr ~ ^^^kr^^kr’^kr^'^kr^^^^^kr^ 

^kr “ ^f^kr^-^kr’-^kr^^kr^'^^^^kr^ 

Then Eq. (2-51) reduces to 


(2-51) 


(2-52) 


E[S|,^(ykr.ykr)^jl = + ^krykr^^j’ 

E[S|^^(;kr.ykr5^j^,.= Et^ykr^kr ^ ^krykr^^j] 


(2-53) 
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Hence, there exists a linear system with spring constant and 
damping co-efficient defined by Eq. (2-52) such that if the 
nonlinear system is replaced by this linear system, the expectation 
values E[g(x,7)x‘] and E[g(x,)0x ] will not be changed. 

Substituting Eq. (2-53) in Eq. (2-50) gives 


r 

r.+=k 

E[g,<(?.3r)Xjl = ecz(y|,^K^ + 
r=fk 


(2-54) 


Let the stiffness matrix and the damping matrix of the linear system 
defined by Eq. (2-52) be denoted by and Then Eq. (2-54) 

can also be written as 


n (e) (e) 

E[g^(x,?)Xj] = E[^E^(C|^3 Xj + Xj)Xj) 

— ' ^ f 0 ) f 0 ) * 

E[ 9 k(x,x)Xj] = E[ £^{C|^5 X5 + K|^5 X5)Xj] 


(2-55) 


since the right-hand side of (2-54) and (2-55) are just two different 
representations of the total force acting on the k th mass. In 
matrix form Eq. (2-55) becomes 

E[g(x,)0x^] = c(e)E[xx^] + K^®^E[xx^] 

_ j_ r \ .:._T / \ 1-T (2-56) 

E[g(x,)0x'] = C^®^E[xxl] + K^®)E[xx'] 

Substituting Eq. (2-56) in Eq. (2-49) and solving for K and C which 

minimizes E[e^e"], yields 

(K - k(®^)E[xx''’'] + (C - c(®Me[xx^] = 0 

(K - k(®))E[xxT] + (C - C^®^)E[xx*'] = 0 


(2-57) 
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which can also be written as 
EfxxT] EExx*"] 

E[xx^l E[xx^] 

If the square matrix is non-singular, the only solution to Eq. (2-58) 
is 


(K - 
(C - 


= 0 


(2-58) 


K = 

C = • 


(2-59) 


If the square matrix is singular, (2-59) is not the only solution. 

But in this case any solution to Eq. (2-58) will lead to the same 
minimum for E[ii^]. So we still can use Eq. (2-59). 

Thus, it has been shown that the linear system formed by 
replacing each nonlinear element by a linear spring and a linear 
damper defined by Eq. (2-52) will minimize E[ee^] provided that the 
excitation is Gaussian. 

It should be pointed out that the smallness of the nonlinear 
term in the original equation was not a presupposition in the develop- 
ment of the above method. In general, the accuracy of this method 
depends on the smallness of the nonlinearity. When the nonlinear 
terms are not small, this method is often an iterative procedure. 

This concludes the relevant literature survey. In the light 
of the review study in this chapter, we will be better able to 
appreciate, in the next chapters, the power and usefulness of the 
simulation approach in the analysis of random linear and nonlinear 
problems . 



CHAPTER III 


DISCUSSION OF SIMULATION METHODS 


In this chapter we will present two methods to simulate a set 
of Gaussian processes {f-|(t), f 2 (t), ... ••• fM(t)} when the 

power spectral density function of each process and cross spectral 
density functions between any two processes are known. These 
spectral densities are usually arranged in the form of a matrix 
commonly known as spectral matrix. This matrix is given by 
Fg^iCw) G]2(w) ... Gin|(w)l 


G(w) 


G2i(w) G22(w) ... G2h/i(w) 


(3-1) 


I^G^^l (w) G|^ 2 (w) ... G^m(w)J 

where G^p(w) is the one-sided cross spectral density between the pro- 
cess fm(t) and fpCt). When m +n, (^n(w) is complex with real part 
Cmn(w) called the co-spectrum and imaginary part - Qmn(w) called the 
quad^spectrum. When m = n, the quad-spectrum is zero and then 
stands for power spectrum for that particular process. 

I . Trigonometric Model 

Here we will show that a set of random processes can be simu- 
lated by a trigonometric series of the form 


N 

f^(t) = Z {a-|^(k)cos [W|^t + aTi(k)] 
k=l 

+ b^^(k)sin iWj^t + a^-|(k)]} 


(3-2a) 


22 
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N 

f2(t) = Z {a2](k)cos [Wj^t + ot2i(k)] + b2i(k)sin [W|^t + a2i(k)] 
k=l 

+ {a22(k)cos [W|^t + a22Ck)] + b22(k)sin [W|^t + a22(k)]> 

(3-2b) 


N 

fj^(t) = z {a^^(k)cos [W|^t + Oj^iCk)] + bj^^(k)sin [W|^t + 0^,1 (k)] 
k” 1 

+ afj, 2 (k)cos [W|^t + 0^2^'^^^ bjjj 2 (k)sin [w,^t + 

+ . • . 

+ + o^(k)] + b^(k)sin [w^t + o^(k)] . (3-2c) 

More generally, the above set of processes can be represented as 
m N 

fm(t) = 2 E {aj^p(k)cos [Wj^t + Oj„p(k)] 

P** I K*“ I 

+ bf„p(k)sin [W|^t + Oj„p(k)]. (3-3) 

The coefficients a^^pCk) and bj^p(k) are Gaussian random variables with 
zero mean that have the following additional properties: (i) afj,p(k) 

and bpq(JL) are always independent and for the case when m = n, p = q 
and k = 1 they are identically distributed, (ii) a^p(k) and apq(Ji) 

(or t>j^p(k) and bpq(Jl)) are statistically correlated if p = q and 
k = 5,, otherwise they are independent. 

The phase angles oij^p(k) are deterministic. Since the co- 
efficients a^p(k) and bpq(k) are Gaussian distributed it follows from 
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the Central Limit Theorem that the time series fn,(t) will also be 
Gaussian distributed. 

The form of the simulated time series given by Eq. (3-3) lends 
itself to a very simple physical interpretation. For the first time 
series (m = 1 ) the index of summation p takes the value of unity. 

For the second time series (m = 2) the index p takes the values 
unity and two. The value of p = 2 adds to the second series two 
terms that are independent of the first series. Physically it adds 
some energy to the second process independent of the first process. 
Likewise for each new process two additional terms are tacked on which 
allow that process to have some energy that is independent of all the 
previous processes. 

In order to see how the co-efficients a„,p(k) and bp^Cfi,), and 
the phase angles c%p(k) should be detemined it is necessary to con- 
sider the cross correlation between the time series ffn(t) and fpCt). 
Here for convenience we will assume m £ n. The cross correlation 
between these time series is given by 
“ E[f,^(t)fp(t +t)] 

m n N N 

= E[ Z E Z Z 

p=l q=l k=l Jl=l 

{a^p(k)a^q(a)cos [W|^t + o^p(k)]cos [Wj^(t + t) + 

+ anip(k)bpq(«,)cos [W|,t + On,p(k)]sin [Wj^^(t + t) + apq(Jl)] 

b^p(k)apq(A)sin [W|^t + c^np(k)]cos [w^(t + x) + apq(£)] 

+ bn,p(k)bpq(A)sin [W|^t + On,p(k)]sin [w^(t + x) + apq(£)]}] 

(3-4) 
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The expectation operator can be brought inside the summations, 
and since the expectation of a sum is the sum of the expectations, we 
get 


Rmn(T^) 


m n N N 

Z E Z Z 

p=1 q=1 k=l £=1 


{E[a^p(k)a^q(5,)]cos [W|^t + otn,p(k)]cos [Wk(t + x) + anq(i-)] 

+ E[an,p(k)bnq(5-)]cos [W|^t + oin,p(k)]sin [Wj^(t + x) + ctnq(^)l 

+ E[b^p(k)anq(il)]sin [W|^t + On,p(k)]cos [w^^(t + x) + apq(Jl)] 

+ E[b„,p(k)b„q(S,)]sin [W|^t + Ojnp(k)]sin [w^(t + x) + a^qC^)]} 

(3-5) 


Because amp(k) and bpq(£) are always independent the middle terms in 

the above expression are both zero. Furthermore, since a^pCk) and 

3nq(^-) independent unless p = q and k = Jl, 

Eq. (3-5) reduces to Rmn(^) 

m N 
Rmn(T) = Z Z 
p=l k=l 

{E[a„,p(k)anp(k)]cos [vi^t + a^p(k)]cos [W|^(t + x) + cipp(k)] 

+ E[bmp(k)bnp(k)]sin [w^t + oimp(k)]sin [wk(t + x) + anp(k)]} 

(3-6) 


^^^’^mnp^’^^ "" E[amp(k)anp(k)]. That isK^np(k) is the covariance 
between a„,p(k) and app(k). Because bj„p(k) has the same distribution 
as amp(k), we also haveK^p(k) = E[bj,,p(k)bpp(k)] . Using these 
relationships and a trigonometric identity Eq. (3-6) becomes 
m N 

Rfrin(^) - ^ ^ “ ®*mp(k) oinp(k)] 

p=l k=l 


(3-7) 
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The right hand side of Eq. (3-7) is independent of the time t; there- 
fore, the assumed trigonometric form for the time series yields 
stationary random processes. 

Noting the fact that the spectral density is the Fourier Trans- 
form of the correlation function, the one-sided cross-spectral 
density for processes fni(t) and fp(t) is given by 

r 


^mn(w) - 


r _ , ,-jWT 

R^„(t) 6 dx, w > 0 


(3-8) 


1^0, w < 0 

The next step is to substitute Eq. (3-7) into Eq. (3-8) and perform 
the integration. However, it should be recalled that Kmnp(l^) Eq. 
(3-7) is the expected value of a^p(k) times a^p(k) (or b^p(k) times 
l3np(l^))* Therefore, by substituting Eq. (3-7) into Eq. (3-8) what 
we obtain is an expected value for the cross-spectral density. Thus 
we have 


E[G|nn(w)] = 


r 


m 

2 S 
p=l 


N 

S 

k=l 




“jWT 

cos [W|^T - ctfnp(k) + anp(k)]e dx} 


w ^ 0 

0, w < 0 

Carrying out the integration gives 
m N 


(3-9) 


E[Gmn(w)] - ] 


2tt 


)=1 k=l 


X {cos [oin,p(k) - anp(k)] - j sin[on^p(k) - app(k)]} ,w > 0 
^0, w < 0 (3-10) 

where 6(w - w^) is the Dirac delta function. 
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We thus see that for the assumed trigonometric time series given by 
Eq. (3-3) is composed of both real and imaginary terms. When 
m = n, Gnin corresponds to a power spectral density and the imaginary 
portion is zero. A schematic plot of Eq. (3-10) is shown in Figure 
1. The frequencies W|^ are chosen such that 

W|^ = Wj^ + (k - ^)Aw, k = l,2,...N 

AW = V (3_n) 

N 

where w^ and w^ are respectively the upper and lower cut off fre- 
quencies of the spectra that is used to generate the time series. 
Suppose an actual power spectrum 6(w) for which we want to obtain a 
simulated time series is shown in Figure 2. We can slice this 
spectrum up into N intervals of width Aw as shown. The amount of 
energy contained in the slice centered at w = W|< is G(w|^)Aw. We want 
the expected value of the spectrum corresponding to our simulated 
time series to have the same amount of energy at w = wj^. This means 
that the expected amplitude of the delta function at w = w|< should be 
G(w)Aw. Similar criteria must be applied for the cross-spectra of the 
simulated time series. 

We assume that at any frequency w = Wk the cross -spectral 
matrix can be factored into an upper and lower triangular matrix, 
where the upper triangular matrix is the complex transpose of the 
lower triangular matrix, as indicated below: 
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Fig. 2. Schematic sketch of a continuous power spectrum 
used to generate a time series 
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Gii(W|^) 

6i2(W|<) •• 

. .. Sif^(W|^) 

G 21 (W[^) 

S2^'^k^ •* 

••• 

Si 

6(^2 

... 


H.|.|(W|^) 0 ... 0 


Hii(W|^) H2 i(W|^) ... Hmi(W|^) 

H2i(Wk) ^22^S^ *** ^ 


0 ^22^^^k^ ^M2^'^k^ 

Hf^l(Wk) S2^'^k^ 


0 0 


(3-12) 


The bar is used to denote the complex conjugate. 
If Eq. (3-12) holds then 

n _ 

Gmn(w) = 2 H (w)H (w) 
p=l 


(3-13) 


for n <_ m. Eq. (3-13) can be rewritten with the limit of summation 
taken as m instead of n since H^p(w) = 0 for m > n. Thus Eq. (3-13) 
becomes 


m 


G|iin(w) " E ^mp^''^^^np^''^^ 

P=1 

Using the above relation,solutions are obtained as 


(3-14) 


CL(Wk) 


H. 


, m = 1 , 2, M (3-15) 

where I^(w(^) is the m th principal minor of G(w|^) with Dq being 
defined as unity, and 
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Hnp(wk) = HppK) t ::::: g : ? 



Dk(w) 


(3-16) 


P = 1 . 2 

M, m = p + 1 , . . . 

. M 

where. 

Gii • • 

. . • • Gi ,p_i 

^Ip 

G/ 1 , 2, . . . , p-1 , m\ — 
M, 2, ..., p-1, p' 

G 21 G 22 • • 

• * • • *^2,p-l 



Vl ,1 . Vl .2 

• • • • S-1 »P-1 

S-1 .P 


^1 * ^2 • * 

• • * • ^,p-l 

^p 


is the determinant of a submatrix obtained by deleting all elements 
except (1, 2, ... p-1 , m) th rows and (1, 2, ... p-1 , p) th columns of 
6(w|<). 

It is noted that the above solutions are valid only when the matrix 
G(w) is Hermitian and positive definite, as can be seen from Eq. 
(3-15). 

Because the cross -spectral density matrix G(w) is known to be only 
non-negative definite, special consideration is needed in those cases 
where G(w) has a zero principal minor. This is discussed in Appendix 

C. 

For the band of frequency [W|^ - ^ Aw, W|^ + ^w] we want the ex- 
pected value of the energy in the spectrum of the simulated time 
series to be equal to the energy Gmn(w)Aw of the original spectrum. 
Thus equating the energy in the band centered at w = W|^ as given by 
Eqs. (3-10) and (3-14) we get 

2TTK^np('<)^cos [On,p(k) - Onp(k)] - j sin [a^p(k) - oipp(k)]} 

= Hmp(wk)H„p(Wk)Aw 


(3-17) 
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If we write 

= |Hmp(w|<)|[cos a^p(k) - j sin On,p(k)] (3-18b) 

= RE[Hn,p(w,^)] + j IM[H^p(Wk)l (3-18c) 

then by substituting Eq. (3-1 8a) into Eq. (3-17) we can see that 

“ ^l*^mp^''''k) 1 l^np^'^k^ I (3-19) 

From Eqs. (3-1 8b) and (3-1 8c) we find 


%p(k) 




IM[Hmp(wk)] 

RETHjnpTwj^TT 


(3-20) 


Thus it follows, as stated earlier, that the phase angles are deter- 
ministic. 

We recall that 




E[amp(k)a 

nn (k)] 


np' 


or 

E[b^p(k)bj^p(k)] 


(3-21a) 

(3-21b) 


That is we want a,np(*^) independent random variables 

with their covariances given by the right hand side of Eq. (3-19). 
This can be done in the following manner, if we let 


_ Aw I 

2Tr'"nip' 


^mp^'^^ " PTfl^mp^'^k^ 


and 




(3-22a) 


(3-22b) 


M and E and n are in- 
P P 


f"P' ' 2 tt' ™p'"K''”p 

where m=l,2, M;p=m,m + 1, 

dependent Gaussian random variables generated on a digital computer. 
Because the ^'s and n's are Gaussian distributed, the a's and b's 
will also be Gaussian distributed. 
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Substituting Eq. (3-22a) in Eq. (3-21a) or Eq. (3-22b) in Eq. (3-21b) 
and taking the expectation of both sides we see that in order to 
make the left-hand-side equal to right-hand-side, the mean of ^ or 
n should be zero and their variances should be unity. 

From Eq. (3-21 a) we have 

E[a„p(k)l = ^|H„p(wk)iE[EpJ = 0 (3-23a) 

and 


" ^I^mp^''^k) 1 l^np^^'^k^ I ^t^p ^ 

(3-23b) 


Thus the terms a„,p(k) (or b^p(k)) satisfy the properties stated earlier 
and they also satisfy Eq. (3-17). 

This concludes our explanation of how to simulate several multi - 
correlated time series. 

After discussing the Fast Fourier Transform in the next section, we 
will give an example problem to check and compare the two methods. 


II . Fast Fourier Transform (FFT) Method 
Assume that a set of discrete time series can be represented by the 
following expression: 




N-1 
E X 
k=0 


Pk 


EXP 


2Trkn 

N 


(3-24) 


where , 

p = number of time series 
n = time index 
k = frequency index 
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The terms Xp|^ are elements of a complex vector generated in the 
following way 

where. 


h = At, discrete time increment 

[H|^] is a lower triangular matrix obtained from the given 
spectrum matrix as discussed in the preceding section, that is, it 
satisfies the relation 

[HklEH^]'^ = [Gkl (3-26) 
and, 4p|^ - ^p[^ + j ripj^ (3-27) 
where and rinu are independent Gaussian random numbers such that 

(3-28a) 
(3-28b) 


E[?pkl = Efnpk] = 0 
ElSpkl = E[n2|,l = 0.5 
Writing Eq. (3-25) in expanded form 


^Ik 

^2k 

i. 

Hllk ° .... 0 

^21 k ^22k • • • • 0 


^Ik 

C2k 


fN f 
" I2h 




^Mk 


^Mlk ^M2k • • • ^ HjviMk 


CMk 

L J 


That is, we can write 




= [zh] 


Vj- 


N C 

2hjq^l”pqk^qk 

N_1^M _ _ 


where the bar denotes complex conjugate. 


(3-29a) 


(3-29b) 

(3-29c) 
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Then, 


M M 


E[XpkX 2h ^ ^pqk*^rsS-^qk^s£^ 

q=l s=l 

Moving the expectation operator inside, we can write 


(3-30a) 


N 


M M 


^f^pkVjl^ 2h ^^i^pqk^rs5,^*-^qk^sil^ 


(3-30b) 


Making use of properties of complex random number as given in terms 
of ^'s and n's by Eqs. (3-28a) and (3-28b), it can be easily seen 
that 

E[?qk?sil] = (3-30C) 

where 6qg and 5kj^ are Kronecker delta functions; substituting (3-30c) 
into (3-30b), we get 


M M 


E[XpkXj,^] = 21^ ^^^*^pqk^rsJl^qs^kil 


(3-31a) 


■ 2h ?i^Pqk^rq£‘^ 


k£ 


Taking the inverse Fourier Transform of Eq. (3-24), we get 

Xpk = fpnEXPI-j 
n=0 


(3-31b) 


(3-32a) 


Therefore , 


Y -Vf FYPr 27r(N-k)n , 

^p(N-k) " ^ ^pn^^*^f 3 |\| ^ 

n=0 


2irkn 

= 2 fpnEXP[j 2^] 

n=0 ^ " 


(3-32b) 
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Since fp^ is real, comparing Eqs. (3-32a) and (3-32b), we see 


that 


(3-33) 


N 


^p(N-k) “ ^pk 

Thus, it is only necessary to generate ^ values of X for each time 
series because the other half of the values of X are just the complex 
conjugate of the first half of the values of X. 

Now, it is necessary to show that time series represented by 
Eq. (3-24) do indeed have the proper power and cross-spectral 
densities . 

From Eq. (3-24) 

N-1 




Since fp^^ is real, taking the complex conjugate on both sides of 
the above equation, we can write 


1 27tnk 

fpn 


k=0 


N 


(3- 34a) 


Thus, the cross correlation between time series p and time series r 
is given by 

Rp^(n,m) = E[fpnf^l (3-34a) 

where, n and m stand for the time indices of process p and process r 
respectively. 

Substituting for fp^ and f^ from Eq. (3-24) and Eq. (3-24a) 
respectively in Eq. (3-34a), we get 

Rp,(n,m) = V 


(3-34b) 
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Moving the expectation sign inside 

1 N-1 N-1 

Rp^(n,m) = ~2 

X EXP[j 

From Eq. ( 3-31 b) one obtains 

— _ li ^ — 

■ 2h q^^^pqk^rqJl'^kfi,. 

Substituting this in Eq. (3-34c) 

1 N-1 N-1 M 

Rpr(n,m) = ^ ^^^^^pqk^rqJl'^kil^ 

X EXP[j (3-34d) 

Unless k = )l, Eq. (3-34d) is zero. 

Hence, 

Rp,(n.m) = ?I^1. (3-34e) 

That is, the correlation is a function of the time lag which 
proves that the process is stationary. 

Finally we can write Eq. (3-34e) as 

= m ^2^,. (3-36) 

The spectral density function is the inverse Fourier Transform of 
Eq. (3-35). That is 

Gprk = 2h VRpr(m-n)EXP[-o MfelL] (3-36a) 

M 

^/^pqk^rqk 
q=l 


(3-36b) 
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which is equal to the elements of the spectral matrix between the 
process p and process q. 

Thus, we see that the assumed form of the time series given by 
Eq. (3-24) has the same spectral density as the target spectral 
density function. 

Simulation of Strong Wind Turbulence 

In order to check the two simulation techniques described above 
namely the trigonometric model and Fast Fourier Transform model, we 
generated six correlated time series representing the fluctuating 
part of wind velocities. These six series represent the three com- 
ponents of the wind at two different spatial locations. Each simu- 
lation technique was checked by comparing the spectra of the simulated 
data to the spectra used to generate the data. The input spectra 
used for these time series are based on wind velocity measurements 
made by several investigators [16, 17, 18]. These spectra will not 
be discussed here in detail. Under strong wind conditions turbulence 
is due to mechanical mixing caused by frictional effects. In this 
case, the turbulence can be considered a stationary random process with 
zero mean and a Gaussian distribution. 

The following expressions for the spectra were used for this 
simulation: 

nC-|-|(n) 87.56f-| 

u*2 ({1.0 + 1.5 (19.6fi)°'®^^)l 

nC22(n) _ 24.23fi 

u*2 ((1.0 + 1.5 (7.368fi)^-^®^)2-134 


(3-38) 
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nC33(n) 

u*2 

3.36fi 
1 + lOfi^^^ 


(3-39) 

nC44(n) 

272. 2f2 


(3-40) 

u*2 

\ - U.S45 1 

((1.0 + 1.5 (39.3f2) ) 

.972 

nCggCn) _ 

46.46f2 


(3-41) 

u*2 

((1.0 + 1.5 (11.06f2)°-^®^)2*'' 

nCgeCn) _ 

3.36f2 


(3-42) 

u*2 

1 + 10f2^/^ 


nC^3(n) 

12.5fi 


(3-43) 

u*2 

(1.0 + 7.5f])®/2 


nC4g(n) 

12.5f2 


(3-44) 

u*2 

(1.0 + 7.5f2)^/2 


C-i 4 (n) = (Cy|(n)C44(n))^(exp(-19Af))^cos 

(2iTAf) 

(3-45) 

Ql 4 (n) = (Cii(n)C44(n))^(exp(-19Af))^sin 

(2TrAf) 

(3-46) 

C25(n) = (C22(n)C55(n))^(exp(-13Af))^cos 

(47TAf) 

(3-47) 

Q25(n) = (C22(n)C55(n))^(exp(-13Af))^sin 

(47TAf) 

(3-48) 

C36(n) = (C33(n)C55(n))^(exp(-13Af))^cos 

(4iTAf) 

(3-49) 

Q3e(n) = (C33(n)Cg5(n))^(exp(-13Af))^sin 

(47TAf) 

(3-50) 

C34(n) = (Ci3(n)C46(n))^(exp{-16Af))^cos 

(STTAf) 

(3-51) 

^16 " ^34 



(3-52) 

Q34(n) = (C-|3(n)C4g(n))^(exp(-16Af))^sin 

(SirAf) 

(3-53) 

Qi6(i) = Q34(n) 


(3-54) 


where and Q^j represent the co-spectrum and quadrature-spectrum 
parts of the spectral density function respectively. 



40 


Note that C^-j = Cj^ 
and q,j =-Qj, 

The elements of the components of the spectral matrix, which have not . 
been listed above, have been assumed to be zero. The following actual 
values were used for simulation 
z-j = 100 ft, Z2 = 50 ft 
V| = 60 ft/sec., V 2 = 50.974 ft/sec. 
f-| = nz-|/v-|, f 2 = nz 2 /v 2 
u*2= (5.212)2 = 27.1649 ft^/sec^ 
f = nAz/v = n (100 - 50)/(60.0 + 50.974)/2 
where 7 = average velocity 

n = frequency in cycles per sec. 

Suffix 1 or 2 refers to the station at 100 ft or 50 ft respectively. 
Results 

Fig. 3-a shows a set of six simulated correlated time series 
based onthe trigonometric model. Time series one, two and three repre- 
sent the wind velocities in the mean wind, cross wind and vertical 
direction respectively at a height of 100 ft. Time series four, 
five and six are similar simulations for a spatial position 50 ft. 
below the first point. Series one and four are positively correlated, 
while series three and six are negatively correlated to one and four. 
Series two and five, the series representing the cross wind, are 
positively correlated to each other but uncorrelated to the other 
four series. For these series, velocities were calculated every second 
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for 2048 seconds. The number of frequencies used, equal to N in 
the analysis, was 600. 

Fig. 3-b shows a similar set of simulated time series based on 
the Fast Fourier Transform method. It is interesting to note the 
similarity of the outcome. The number of frequencies used, equal 
to N in the analysis, has to be the same as the number of time points 
that is 2048. 

A comparison between the spectra used to generate the simulated 
data and the spectra of the simulated data is presented in Figures 
4-a through 9-a for the trigonometric model and in Figures 4-b 
through 9-b for the FFT model. Figures 8-a and 9-a are respectively 
the CO- and quad-spectra between time series one and time series four 
based on trigonometric model. Figures 8-b and 9-b are the FFT counter 
part. In each case the solid line represents the original spectrum 
and the dashed line represents the spectrum of the simulated time 
series. The spectra were calculated from the simulated data using 
the "Cool ey-Tu key" method of calculating the correlation function 
first, and then the "Tukey window" was used for smoothing purposes. 

As can be seen the comparison between the actual spectra and the 
simulated spectra is quite good in all cases. We did not consider it 
necessary to show all 36 possible power and cross-spectra. The 
agreement looks the worst for the quad-spectrum between series one 
and four (Figures 9-a and 9-b). The vertical scale for this compari- 
son, however, is very expanded and about all we can conclude is that 
both the original and simulated spectra are very close to zero. Also, 
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Fig. 6b. Comparison between the original G 33 (w) and the 

spectrum of the third simulated time series based 
on FFT model. (Vert. Scale: 5(M2/sec)/in; Horz. 
Scale: 0.63 6"ad/sec,)/in) 
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Fig. 7a. Comparison between the original G 44 (w) and the 

spectrum of the fourth simulated time series based 
on trigonometric model . (Vert. Scale: (150 M^/sec)/ 
in; Horz. Scale: 0.63 (rad/sec.) /in) 
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Fig. 7b. Comparison between the original G 44 (w) and the 

spectrum of the fourth simulated time series based 
on FFT model. (Vert. Scale: 150(M^/sec)/in; Horz. 
Scale: 0.63 (rad/sec )/in) 



2H'3 32-» 99 


Comparison between the original C-|^{w) and the co- 
spectrum between simulated time series one and four 
based on trigonometric model. (Vert. Scale: 

81 .25(M^/sec)/in; Horz. Scale: 0.63. (rad/sec)/in) 
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Fig. 8b. Comparison between the original G-]^(w) and the co- 
spectrum between simulated time series one and four 
based on FFT model. (Vert. Scale: 81 .25(M^/sec)/in; 
Horz. Scale: 0.63 (rad/sec )/in) 
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Fig. 9b. Comparison between the original Qi 4 (w) and the quad- 
spectrum between simulated time series one and four 
based on FFT model. (Vert. Scale; 4.0(M2/sec)/in; 
Horz. Scale: 0.63 rad/sec) /in) 
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it is interesting to note that, in general, the comparison between 
the simulated spectra and the actual spectra looks better in the 
case of FFT approach. 

Comment 

Considering the above results we feel that both the simulation 
methods presented work very well. The only drawback we see is the 
amount of computer storage that is necessary. It is noted that FFT 
method works about twenty times faster in terms of computer time 
compared to trigonometric model. But FFT approach has a drawback in 
that the number of time points is equal to the frequency points and 
we can't restart the simulation at any given time t. It has to 
start with t =0. Also, the time interval is related to frequency 
intervals. The trigonometric approach does not suffer from these 
disadvantages. Here, we can start the simulation beyond any given 
time and also the time interval is not dependent on the frequency 
discretization. Because of the presence of ^ine and cosine functions 
in the trigonometric model, it can be very useful in the response 
analysis of linear system. Thus, we see that whereas the Fast 
Fourier Simulation method has a tremendous advantage in terms of 
computer time, trigonometric model is more useful in certain cases. 



CHAPTER IV 


APPLICATIONS OF THE SIMULATION TECHNIQUE 

Example 1 . Forced oscillations of a free standing latticed tower. 

To show the application of the simulation method for a linear 
problem the dynamic response of a three-legged, latticed steel 
tower, shown in Figure 10, due to wind load in the y direction was 
investigated. The tower was idealized as a space truss, and the 
structural members were assumed to resist axial loads only. 

For the dynamic analysis, the tower was modeled mathematically 
as a discrete system of fourteen masses lumped at the panel points. 
Horizontal loads were assumed to be applied at the panel points, and 
secondary stresses were assumed to be small. 

The assumptions made in deriving this analytical model are 
summarized below: 

1. The tower is a linear elastic space truss; 

2. Motion in two orthogonal horizontal directions are uncoupled; 

3. Vertical motions are negligible; 

4. Loads are applied only at panel points; 

5. Secondary stresses are negligible; 

6. Masses are concentrated at centroids of planes at panel 
points; 

7. The base of the tower is rigid; and 
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8. The deformations are small enough to be geometrically linear. 
Under these assumptions, the column elements of the flexibility matrix 
were computed by successively applying a unit load horizontally at 
each panel point in y direction. The flexibility coefficients are 
shown in Table 2. The mass matrix is diagonal . [37] 

After the flexibility and mass matrices for the tower had been 
computed, the natural frequencies and mode shapes were obtained in 
the usual manner as described in the following section. 

The equation of motion for an undamped, multi -degree of freedom, 
discrete mass systems are given by 

My + Ky = F(t) (4-1) 

where M = mass matrix, 

K = stiffness matrix of the system, 

F(t) = forcing function vector of the system, 
y = displacement vector, 

and, 

y = acceleration vector of the system. 

When the system is vibrating harmonically at a natural frequency, 
F(t) is zero and the displacement vector may be written as 

y = ApSin(wnt) (4-2) 

where A^ = displacement vector of the n-th mode 
Wp = n-th natural frequency. 

Eq. (4-1) then becomes 

-Wp^MApSin (Wpt) + KApSin (wpt) = S’ 


(4-3) 
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2 _.. 


Dividing Eq. (4-3) by Wp sin (Wpt) and rearranging terms, then 


MAn = 


w, 


n 


(4-4) 


Premultiplying both sides of Eq. (4-4) by the inverse of K 
K"^MAp = K“^KAn 

Wn^ 

and replacing K~^ by S, the flexibility matrix of the system, and 


(4-5) 


K~^K = I = identity matrix, then 
SMAp = Ap 

Wn 


(4-6) 


1 


Denoting D = SM as the dynamical matrix and Ap = — standard 

Wp 

algebraic eigen value problem is obtained as 

D^p = ^n^n (4-7) 

The computed first three eigen values and their corresponding eigen 
vectors are shown in Table 3. Three natural periods were 0.496 sec., 
0.158 sec., and 0.080 sec. in the y direction. 

After having determined the modal shape and natural frequencies, 
we will proceed to do the dynamical analysis as follows: Let the 

masses Mj, concentrated at points 1, 2, 3, ... j, k, ... r of the 
system be subject to disturbing forces. 

(4-8) 

where = O-SpCxjSjVgj -- the static wind loading acting at the 
point under consideration, v^j is the average value of the longitudinal 
component of the velocity at this point, 

Foj(^) ~ ^ ) is the disturbing force relevant 


= Fgj + F^j(t) 


aj 


'aj 



for the pulsating part of the velocity VQj(t), 
p is the air density, 

is the drag coefficient of the construction at the level j, 
and Sj is the area of the projection of the structure at the level j 
on the plane perpendicular to the direction of the wind. 

Now, let 


1=1 


(4-9) 


where y|^ = dynamic displacement 

c^^-(x|^) = modes of natural oscillation 
q^-(t) = generalized coordinates. 

Then the equation of motion for the i th mode is 

q.(t) + (u + iv)w^-^q^(t) = SlilL (4-10) 

Hi g 


where u = 1 


ro 


IV = 


4 + ro^ ’ 
4rp 

4 + rn^ 




The generalized force Q-j(t) = 


2 

mass M-jg = I l^'Ij'^yi ^^j^ where 
J 1 


r 

E Fj(t)ot -j (xj) , and the generalized 
j=l ^ 


Wj = i th natural circular frequency of the system 
6 = logarithmic damping decrement of oscillation 
Neglecting the square term in the expression for the disturbing 
force, which is very small compared to the first term, we can write 
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Now, 


0 - 


'aj 


Then 


Q^-(t) = mean generalized force = EFajotyj(xj) 
Bq-j(T) = correlative function of the generalized 
force = E[Qj(t)Q^-(t + t)] 


N N 


= E[ E E {F (xJ(1 + 

j=l k=i aj yi J Vgj 


(4-12) 


(4-13) 


X ^Fakayi(xk)(l + (4-i4a) 

Va 1 


'aj 

Multiplying these two expressions and moving the expectation sign 
inside 


N N 
N N 

^ yi(^k^ 


^ E[ Voj(t)Vok(t + t)] 

''aj^ak 

Bq^qj^(t) = Cross correlative function for the generalized force = 
E[Qi(t)Q^(t + t)1 

N N 2v -ft) 

X {F,kc.yj(xk)0 (4-15a) 
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N . 

®QiQ5, ■ ^ 

J ” I K I 


N 

+ Z 

j=l 


N 

1 ^^ 1 j ^ak“yi ^ ^ j ^ “yJl ( ^k ) ^ [ 


Vi(t)Vok(t + t) 
''aj^ak 


■] (4-1 5b) 


Taking the Fourier Transform of the correlative and cross correlative 
functions we get 

Sq^(w) = spectrum of the i th generalized force 
N N 
j=l k=l 

f° (1 + 4 ^^^oj^^)Vok(t ■» t)] p-iwT^^ (4-16a) 

'^aj^'ak 

Assuming the mean wind to be invarient with time 

N N 

SQi(w) = Z 2 FajFg|^ayj(Xj)ay^'(X|^) 

J *" I K"* I 


X (6(w) + 


''aj'^ak 




(4-16b) 


where SVojVgk is the cross-spectrum of the fluctuating velocities at 
points j and k. Note that Sq-|(w) will be complex because SvojVqi^ might 
have real and imaginary parts. 

However, the imaginary part can be neglected if this turns out 
to be very small compared to the real part. 

Similarly, 

SQiQ£(w) = cross -spectrum of the generalized force 
N . N 

^QiQ£(w) = Z ^J^aj^ak<*yi(^j)“yz(^k^ 
j=l k=l •’ 

X <«(“> + SVojVoi;) 


(4-17) 



For considering the dynamic effect, we don't consider the static de- 
flection and in such cases the spectral density for the fluctuating 
part only will be 


N N 

SQi(w) = l^f^^^aj^ak^yi^^j^^yi^^k^ 


and. 




N N 


( 4 - 18 ) 




( 4 - 19 ) 


'aj^ak 

Now, to find: the expression for the transfer function of the system. 


let 

Qi(t) = eJwt 
and, q^-(t) = $j(w)eJ^ 
where $^-(w) is the transfer function. 

Substituting these expressions in Eq. (4-10) we get 

-w2$. (w)e'^'^^ + (u + iv)wi^4>,*(w)e^'^^ = — 

' Mi g 

_ 1 

Mig[-w^ + (u + iv)wi2] 

-w^ + (u - iv)w^-^ 

M^-g(-w^ + (u + iv)w^-^)(-w^ + (u - iv)w^-^ 
-w^ + (u - iv)wj^ 

~ M-jg[w^ - 2 uw^Wi '2 + w-j^] 


or> 4>j(w) 
(w) 

<J>i(w) 


( 4 - 20 ) 
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Since = 1 


M^-g^[w^ - 2 uw^w^- 2 + w^-^] 


(4-21) 


Then, 


Sq-j(w) = spectrum of the generalized coordinate 
Sqi(w) = Sq^-(w)I4>^-(w)|2 

Spi (w) 

Sni(w) = 


' ' M-jg^(w^ - 2uw^w.^ + w^-^) 
and using the relation 


(4-22) 


Sqiqil^w) = Sqi-qj,(w)$,-(w)$^(w) 

we get the cross-spectrum of the generalized coordinate as 
S ^ SQiQ£(w)[w4 - u(w^-^ + Wj^2)w^ + iv(Wj^ - W;^^) ->■ 

M^gMj^g(w^ - 2uw-j2w2 + w-j^)(w^ - 2uWj^^w^ + 


(4-23) 


Now, for the motion of any point k 

r 


y,^(t) = S qi(t)ay|(X|^) 

i=l 


Similarly for any point % 

r 

y^^(t) = E^qs(t)ay5(x^) 

Then correlation function between y|^(t) and yjj^(t) 

r r 

^qiqs^^^*^yi ^^k)“ys^^JJ.^ 

Then the spectral density of the dynamic displacement 


(4-24) 


sykyji = ' 




-oo 
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r r 

S ^ “yi(^k^“ys(^£)5qiqs 
1=1 s=l 

r r _ 

SykYjj, = ^ ayi(X|^)ays(x^)$i(w)$(w)SQ^-Qj;^ 

i=l s=l 

r r 
1=1 s=l 

N N ^SVojVok 

^ l^^^''ajFak“yi^^J^°‘ys^^<^^ VajVak 

[w4 - u(w^^ + Wg^)w^ + iv(w-j2 - Wg^) + 

^ MjgMsg(w^ - 2uw^- V + Wi^)(w^ - 2uWg^w^ + Wg^) 

when k = A, it gives the power spectrum of the displacement. 

Thus, having obtained the expressions for cross and power spectrum, 
we can simulate the displacement at any point. Thus the spectral 
density function of the displacement being known, FFT simulation 
method was used to simulate the displacement of the top of the tower 
shown in Fig. 10 whose important parameters are listed in Tables 1-3. 

The value of the co-efficient of drag was taken to be 0.6 and one per- 
cent of the critical damping was used. 

Fig. 11 shows the simulated samples of the displacement and the 
generalized forces. The first three samples are the time histories of 
the generalized force for the 1st, 2nd and 3rd mode of motion respec- 
tively and the fourth sample shows the time history of the normalized 
displacement of the top of the tower due to all the three modes 
combined together. 


(4-25a) 

(4-25b) 

(4-26) 
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PANEL 

POINT 



SECTION A-A 


Fig. 10. Latticed steel tower subjected to turbulent wind 
load. 



TABLE 1 


CONCENTRATED MASSES AND PROJECTED AREAS 


Panel 

Point 

Elev. 

(ft) 

Projected 
Area (sq. ft.) 

D.L. (lb.) at 
Panel Point 

1 

150 

2.87 

1810 

2 

147 

5.22, 

510 

3 

143 

8.06 

750 

4 

137 

8.75 

810 

5 

130 

10.32 

920 

6 

123 

11.97 

1790 

7 

115 

14.12 

1190 

8 

107 

16.76 

1360 

9 

97 

19.30 

2940 • 

10 

87 

23.32 

1900 

11 

75 

28.08 

3580 

12 

61 

32.35 

2940 

13 

45 

43.45 

5670 

14 

25 

60.37 

5240 


NOTE: 1. Values for projected area are for one face only. 

2. Values for D.L. (Dead Load) are for entire tower. 





TABLE 2 


FLEXIBILITY COEFFICIENTS IN Y DIRECTION 
(inch/kip) 


1 2 3 4 5 6 7 8 9 10 11 12 13 14 


1 0.4113 


2 

0.3890 

0,3702 






3 

0.3599 

0.3441 

0.3231 





4 

0.3190 

0.3066 

0.2900 

0.2652 




5 

0.2762 

0.2665 

0.2537 

0.2344 

0.2119 



6 

0.2382 

0.2307 

0.2205 

0.2053 

0.1875 

0.1697 


7 

0.2001 

0.1943 

0.1864 

0.1747 

0.1610 

0.1473 

0.1317 

8 

0.1667 

0.1622 

0.1562 

0.1471 

0.1366 

0.1261 

0.1140 

9 

0.1308 

0.1276 

0.1233 

0.1168 

0.1092 

0.1016 

0.0930 

10 

0.1006 

0.0983 

0.0953 

0.0907 

0.0854 

0.0800 

0.0739 

11 

0.0708 

0.0694 

0.0674 

0.0645 

0.0612 

0.0578 

0.0539 

12 

0.0440 

0.0432 

0.0422 

0.0406 

0.0388 

0.0370 

0.0349 

13 

0.0223 

0.0220 

0.0216 

0.0209 

0.0202 

0.0194 

0.0186 

14 

0.0062 

0.0062 

0.0061 

0.0060 

0.0060 

0.0059 

0.0058 


Symmetrical 
about diagonal 


0.1020 

0.08430.0736 

0.0678 0.0602 0.0526 

0.0501 0.0453 0.0404 0.0347 

0.0328 0.0302 0.0276 0.0245 0.0209 

0.0177 0.0167 0.0156 0.0143 0.0128 0.0112 

0.0057 0.0056 0.0054 0.0053 0.0051 0.0049 0.0048 
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TABLE 3 

COMPUTED NORMALIZED MODE SHAPES 
(Normalized with respect to bottom panel) 

Y DI RECTI DN 


Panel 

Point 

First 

Mode 

Second 

Mode 

Thi rd 
Mode 

1 

39.93 

-5.60 

1.78 

2 

38.31 

-4.73 

1.23 

3 

36.15 

-3.59 

0.52 

4 

33.00 

-2.00 

-0.41 

5 

29.50 

-0.38 

-1.14 

6 

26.22 

0.93 

-1.56 

7 

22.71 

2.10 

-1.63 

8 

19.47 

2.95 

-1.44 

9 

15.82 

3.61 

-0.90 

10 

12.56 

3.83 

-0.21 

11 

9.18 

3.72 

0.63 

12 

5.96 

3.13 

1.29 

13 

3.19 

2.20 

1.56 

14 

1.00 

1.00 

1.00 



Y(t)(IN.) t ,Q3(t)(KB>S) Q2(t)(KIPS) Q,(t)(KIPS) 



FIG. I I. GENERALIZED FORCES AND DISPLACEMENT 
FOR LATTICED TOWER: 

(A) QCNERAUZCO FORCE FOR FB»T MODE 

(B) QEHERALIZEO FORCE FOR SECOND MODE 
(a QENERAUZEO FORCE FDR THIRD MODE 

(0) NORMAUZEO DYNAMIC 0SF1.ACEMENT OF THE TOP OF THE TOWER 
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Example 2 . Nonlinear Random Vibration of String 

One of the most significant applications of the simulation tech- 
nique is for simulating random generalized forces. The necessity of 
simulating random generalized forces arises when the dynamic response 
analysis is performed in time domain for the purpose of obtaining 
information beyond second-order statistics or when the structure is 
nonlinear and therefore an approximate random response is sought by 
simulating the excitation. In order to assert the validity of the 
preceding discussion, the problem of nonlinear vibration of a string 
will be considered in this example followed by the problem of non- 
linear vibration of a plate in the next example. 

The governing differential equation for nonlinear vibration of 
a string is 


where p is mass per unit length, C is the linear viscous damping, Tp 
is the initial tension, f(x,t) is force per unit length of string, A 
is the cross-sectional area of string, and the boundary conditions are 
u(o,t) = u(L,t) = 0 (4-28) 

Assuming an approximate solution of the form. 


u(x,t) = I bn(t)sin 
n=l *- 


(4-29) 


one obtains 
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M 


Z (pb^ + Cb^)sin 
n=l ' 


= fTo + i5£ n7r.2 . 

2L Jo 

X t-Sb„M„ m . ,mr)2j , 


or. 


(4-30) 


X sin = f(x,t) 

Multiplying both sides of Eq. ( 4 - 3 i) 
obtains 


(4-31) 


by sin !2p. an^j integrating, 


one 


= f^(x,t)sin UHdx 

■'o L 

Dividing by p throughout, we get 

f f(*.t)sfn 


(4 -32a) 


pi 


where 6 = 


(4-32b) 


or, 


2p 


b„ + 26b„ + ( 


'" " . V " 76pl3 

' It / ^(*.t)sin Sp<ix 


(4-32c) 
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Assuming that the fundamental frequency of the corresponding linear 

) 

T 


T “2 

spring is 0.1 Hz or ° = 4 sec , we can write Eq. (4-32c) as 


PL' 


M 


bp + 23bp + [1 + 5 (^)^b|^^]{2mr)^b^ 

^ IC“*1 


8L Z’^r/ 4 .\ ■ nrrrx. 

= ^ f(x,t)sin -r-<lx 

'o in 


(4-33) 

The right-hand side of Eq. (4-33) is the random generalized forces to 
be simulated. 

Let f^(t) = I f(t,X])sin Ill^fLdx^ 

fL 

fp(t + t) = f(t + x,X2)sin ”^^2d x2 
Jo L 


Then , 


Rp,p(x) = correlation function of the generalized force 
Rmn('r) = E[fm(t)fp(t + x)] 

Vf^ f'’E[f(t,xi)f(t + x,X 2 )l 

J n J r\ 


. imrxi , mrX 9 j . 

X sin ___Lsin Adx^dxg 

R„,n(T) = f f R(x,x-| ,X2)sin HlZi^in !!J2^x^dx2 

Jo ■'o L L 


(4-34) 


where R(x,x-| ,X 2 ) is the correlation function of the actual force. 
Taking the Fourier Transform of Eq, (4-34) we get 

S^n(w) = [ [ S(w,x-| ,X 2 )sin !^?3^in (4-35] 

■'o “'o ^ ^ 

where Sp^p(w) and S(w,x-| ,X 2 ) are the spectral density functions of the 
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generalized force and actual force respectively. 
For the particular case we take (Ref. [8]), 

S(w,x-| ,X2) = So(w)e““f'^l 

where 

;.2 Of2 


So(«) = „2^ 


(4-36) 


(4-37) 


a and a are constants and gives the standard deviation of the 
actual force. 

Substituting Eq. (4-36) in Eq. (4-35) 

'o ■'o 

= S„(w)I„„ (4-38) 

'o ■'o 


S||,n(w) = 11 So(w)e *^^sin Hl^iLsin !ip?dxidX 2 


where I„„ = f‘‘ f‘'e-“l«l l*l-*2lsih ll^ln ^^x,dx 2 


or, 


‘mn 


= fsin I^dx- 




p’“l"l*2si„ "™2 


^X2l 

L 


dxi 


= f^in f^XT^2a|w|sin ^alw|(x + T 

L 

fL niTTXi -alwIxi.nTr 
X dx^ + sin -j- — e ' '(L“)dxi 

(aw) + (^)^ 


1 
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1 

(aw) 2 + (J1!L)2 

L 


{alwlL . 6, 


mn 




(aw)2 + (p2 


+ (fllL) (p 


(aw)2 + (p2 

which can be further simplified as 


1 


‘mn 


(aw)2 + (H.)2 


■{a|w|L . 6, 


mn 


/ mrr ^ /mr\ 

+ — ^[1 + {-1)’"'^"+ e’“’'^*‘'((-l)"^^+ 

It 


(«)2 + (p' 


('1-39) 


So, 


.2 Oi 


1 


Uf 

^mn^w^ ^ o - V X — _ P 

a^ + 7T (ctw)^ + (M)2 


{a|w|L . 6| 


mn 


/[mr, /mr X 

+ '' 2 „2 n + (-1)"^" + (-1)"*')]) 

(ow)2 + 

(4-40) 

A numerical example was carried out for the case where 3 = 0.1 x 

2tt sec"^ , a = 4it sec“^ , at = 0.7 sec, 1° = 0.05, L = 25 in. , Tq = 

AE 
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100 lb. , p = ^ Ib/in, E = 3x10^ psi . 
do 

From the derived expression for the spectral density function, the 
generalized forces for 1st, 2nd and 3rd modes were simulated by the FFT 
method. They are shown in Fig. 12. Equations of motion (4-33) for 
n = 1, 2, and 3 were reduced to six first order differential 
equations which were solved numerically by using the Modified Pre- 
dictor-Corrector method. Fig. 13 shows the simulated samples of the 
displacement at the midspan. The first three time series are the 
generalized displacement histories for the 1st, 2nd, and 3rd modes 
of motion respectively. The fourth one gives the actual displacement 
due to all three modes of motion at the center of the string. The 
root mean square aQ of nondimensional response 

i(L,t) = u(t,t)/(t ) 

at the midspan is plotted as a function of the root mean square Of 
of nondimensional forcing function 
f(x,t) = f(x,t)/(^) 

in Fig. 14. The total number of 2000 time points with time increment 
of 0.025 sec. were used for response calculation. For the sake of 
comparison, the response of the linear string was also plotted on the 
Fig. 14. 

Conclusion 

By comparing the results with the one obtained by Shinozuka [8] 
for the same problem, it is found that they are very close. This 
demonstrates the useful Iness of the FFT simulation method. It is very 
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FIG. 12. GENERALIZED FORCES FOR NONLINEAR 
STRING = 

(A) MONDU/ENStONAL GENERALIZED FORCE FOR ICT MODE 

(B) NONOMENSIONAL GEIERALIZED FORCE FOR Zm MODE 

(C) NONDIMENSIONAL GENERALIZED FORCE FOR 3rd MODE 
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N 

1-“ .n 

i‘ 




(B) 


1-0 


•0 

o 






(C) 



NONLINEAR STRING: 

(A) NONDIMENSIOHAL GENERALIZED DISPLACEMENT FOR ItT MODE 

(B) NONDIMEN8IONAL OEMERAUZED DISPLACEIttNT FOR 2 ho MODE 

(C) NONDIMEIWIONAL OEHERAUZEO DISPLACEMENT FOR 3w MODE 

(0) NONDIMENSIOHAL DISPLACEMENT TIME HISTORY AT THE CENTER 

OF THE STRINO 


RMS OF NON-DIMENSIONAL DISPLACEMENT 
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interesting to note that total computer time for 2048 simulated 
points was less than 2 minutes which shows the advantage of time- 
saving due to FFT method. To simulate the same number of points 
using trigonometric model takes about 30 minutes of computer time. 
Example 3 . Nonlinear Panel Response Due to Turbulent Boundary Layer . 

The random vibration of a flexible plate immersed in a fluid 
flow on one side and backed by a fluid filled cavity on the other 
side has been considered in this analysis (Fig. 15). The nonlinear 
plate stiffness, induced in the plate by out-of-plane bending and 
the mutual interaction between the external and internal fluid flow, 
is included. The same problem has been considered by Shinozuka [1] 
where he used multidimensional simulation method to find the deflec- 
tion at the center of the plate. In the present work, the 
FFT simulation approach has been used to do the response analysis at 
the center of plate. The results are compared with the ones obtained 
in Ref. [1]. 

The deflection w of a simply supported plate having geometric 
nonlinearity is described, in dimensionless form, by two partial 
differential equations 

V^w = ^yy + " '^y 

- BW|. - w!^ + ^ ^ ^ (4-41 ) 

= 12(1-v^){w5^j - W^^P (4-42) 


with' boundary conditions 
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Fig. 15. Problem' geometry for nonlinear plate vibration due to 
turbulent boundary layer pressure. 
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w(0,y,t) = w(l ,y,t) = 0 
w(x,0,t) = w(x,l,T) = 0 

= v^j^x,l,t) = 0 


(4-43) 


Wyy ( 0 >y ^ ~ ^ 


where x = 

y = 

w = 
$ = 
V = 

t = 
D = 
E = 
3 = 

P = 


x/a, a = plate length 
y/b, b = plate width 
w/h, h = plate thickness 

nondimensional airy stress function for membrane stresses 
Poisson's ratio 

nondimensional time = ^1/— where 

aby p 

Eh3/12(l-v2) , plate stiffness 
modulus of elasticity 
damping co-efficient 


hD 


p(x,y,T) where p(x,y,t) is turbulent pressure 



p®(x,y,T) where p®(x,y,T) is fluid pressure due to 
external flow 


— a^b^ —Cl 

P(. = ^(x,y,^ where Pc(x,y,t) is fluid pressure due to 

cavi ty f 1 ow 

In solving Eqs. (4-41) and (4-42), the plate deflection w is expanded 

in terms of the normal modes of the corresponding linear plate: 

w^(x>y.t) = T. I b^p(T)sin imr'xsin n7ry (4-44) 

m n 

where b^^(t) is the modal amplitude. 
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After the stress function $ is evaluated [19], the assumed modal 
solution is satisfied in the Galerkin sense by computing its integral 
which is average weighted in turn by each term of Eq. (4-44). The 
result is a system of simultaneous nonlinear integral-differential 
equations for b_(t) which involve the following generalized forces 


_ a2b2 f> (■' _ _ 

p(x,y,t)sin imrxsin mryd^dy 

“'0 ^0 


r=C 


'oo'^oo j Q j Q 

1 ,1 


Finn(^ 


hp«a^“ ■'0 ‘'0 


n p^(x,y,t)sin imrxsin mryd-d- 

n ^ y 


(4-45a) 

(4-45b) 

(4-45c) 


where pto = free stream density 

Uoo = mean external flow velocity 
a^ = velocity of sound, cavity flow 
To promote computation ease, a two mode approximation corresponding , 
to X coordinate and one mode approximation corresponding to y co- 
ordinate to Eq. (4-44) has been assumed. Then the solution for $ [19] 
satisfying Eq. (4-42) is given by 


TT^Eh^ 


$ = 


2“ y{[(l+va2)bii^ + (4+va^)b2^^]y^ 


16D(l-v'^)a' 

+ [(a^ + v)bii^ + (a^ + 4v)b2-| ^]a^x^} 

Eh^a2, t>ll2 ^ _ _ 

+ • { [~g cos 2irx - bn b 2 i cos irx 


+ ^11*^21 


9 

9bi 1 b/ 


cos Srrx + 


b2l2 


32 




cos 4irx] + [ 


8a" 


, --’n"*2i ^11^21 — 

,cos irx - - ^ 2 ^ 2 *^°^ 3irx]cos 2iry} (4-46) 


(l+4a^)^ 
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where ct = ^ 

Substituting Eq. (4-46) and Eq. (4-44) into Eq. (4-41) we get the 
following set of equations 

bi + 3-|b^ + (C-|Q + C-jibi^ + Ci2b2^)t>ii 
= HIT. +>Cf,c + 

■ » • I - • • I ••• I / 

b2 + ^2^2 ^ ^^20 ^ ^21^1^ ^2262^)^21 

= 4(F2 + ACF2^ + XF 2 ®) 

where b-| = b-| -| , b 2 h b 2 i , F-] = F-|^ , ^2 h F 2 i , etc. 

^10 = 

r - / X 4^2 4 

Znn - (a + -) 7T 


r 3 4 
^11 - 4 ^ 


a 


[(l-v 2 )(q 2 + ?) + 2 ( 2 v + + ^)] 

a ot 


C,2 = |it''{0-v2)[4(c. 2 + ij) + 

^ (9m5T2> * 


C21 - C^2 


C 22 = |’>'^[(1-v^)(a^ + 

+ 2(8v + + ^)] 


(a.-in?>'\ 


(4-47b) 


3-j and $2 the damping coefficients for the first and second mode 
respectively. The generalized aerodynamic force of the external flow 
is given by 


C = :: Q 


rm 


where is given in references [20] and [21]. For high mach 


(4-48) 
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numbers, Eq. (2-48) can be approximated to 


e 


r 


(4-49) 


where 


rm 


2M(m^-r^) 
= 0 

r/D 


[1 - (-l)^*"] for r + m 


for r = m 
for r = m 


™ 4MbUx»^ 

= 0 for r=f m, M = mach no. 

Again based on Ref. [21], the generalized cavity force is given by 

(4-50) 


rc_ _a_ rn.U_ :^- l )"! ]. [l_-_(ji1 i :ib 

*^m - - _4^ _,2 I2 




\mi n _ 
m*-’ y.=l r 

Thus having obtained the approximate expressions for the generalized 
forces due to external flow and cavity flow, we proceed to find the 
generalized forces due to turbulent pressure fluctuations as follows 
by usingthe FFT simulation technique: the semi -empirical formula of 

cross -spectral density for subsonic boundary layer turbulence is 
given by [22] 

S(w,k-| ,k 2 ) =-jl 


0.715 X lO"® 3^ 


X [3.7exp(-2j5^) + O.8exp(-O.47j5^) 

Uoo Uoo 


- 3.4exp( 




[(o.lg-)^ + (r + k^)^]l(0.715g-)2 + k/] 


(4-51) 


for 0 < w < <» 

where, w = freq. in radians, k-| and k 2 are wave numbers in the x and y 
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directions respectively. = 0.65 Uoo, q = dynamic pressure = 
2 

pUoo /2, c = boundary layer thickness. 

Let = ^2 - x-| 

?2 = yz - yi 

Thpn +-Vip cr>Q/"hv»al Honci-hw •FimP'Mrm 

VI ^wv»wwi\Ai V4WIIW1 wjr I vtvtwwiwSS 

S(w,C] .?2^ "" f ,k2)dk^dk2 

"00 —00 


TT^lLo ^-°° ' 

X r e'^^2*^2(^L^l§!l)dk2 
i(a^)' * k/] 


where Cq = 3.7exp( 


2 , ,0.1w , ,0»715w _ 

C xl0"^(3^)e"l^ll g 

° tt'^Uoo 

-2^)+ 0.8exp{-0.47^) - 3.4exp(-8^) 


j|wUi 


(4-52) 


Now, the generalized force due to turbulent pressure fluctuation is 
given as 


a2b2 fl fl 

(t) = J J p{x^,y-,,t)sin m-|TTX^ 


mim2 


0 ''0 

sin ni27ry,dx-] ,dy^ 
a 2 fa 2 f"* 


(4-53a) 


— — a^b*^ f' f' — — — _ 

^n.nJ-^ + ^) = ~W P(x 2 .yz*t + T)sin ni7rX2 
• ^ ^ 0^0 


sin n 2 iry 2 dx 2 dy 2 


(4-53b) 


X sin miTTX-jSin 012^1 sin n-iTrx2sin n^iry^dx-j dy-j dx2dy2 


or. 


1 ,1 ,1 ,1 


/„ lo I •^2-^) 

X sin m-iirx^sin m2TryiSin n-jir 72 sin n2iry2dx-] dx 2 dy-| dy2 

(4-54) 

where Rin^m2nin2^^^ cross correlation function of the generalized 

force and Rp()Ti ,X2,y^ ,y2»^ is the cross correlation function of the 
turbulent pressure fluctuations. 

Taking the Fourier Transform of both sides of Eq. (4-54) we get 


mim2nin2 




— - 

S(w,C,,e 2 ) 

0 Jo 


X sin m-iTTXisin m2TryiSin niirx2sin n2iTy2dxidx2dyidy2 

(4-55) 


where S^^ni2nin2^^ S(w,?1,C2) are the spectral density functions 
of the generalized forces and turbulent pressure fluctuations 
respectively and w is the nondimensional frequency. 
Nondimensionalizing Eq. (4-52) with respect to w, C] and ^2 and sub- 
stituting in Eq. (4-55) for S(w,T^ ,f^) we get 
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\2u^^ibVr° 



l^r 




_ ,0.715wb 



W 


0 •'0 •'0 '0 
,aw 


.rt*«»rt/_ e "I n r>_TrV-. 
iii^iijr ^ o I II **^'*^2 


„ g-3(x2-xi)TF... 

A <= ■'a I II nil iiA-j a I II 

X sin n2iry2d7] d^2^yi ^y 2 (4-56) 

In order to solve the integral in Eq. (4-56) we proceed as follows: 
let 

, , 1 — — lO.lwa . — ,aw 

^ -Ix2-^ll“13r -j(^2-xi)lF 


'"■l"! ' lo i 


^2^2 


then 


%in-| 


sin miTTXiSin ni7rx2dxidx2 

,1 rl 

e’ |y2-yi I DT^ 

0^0 

sin m 2 iTy^sin n 2 iTy 2 dy ■] dy 2 

J’UF^I . -O-rjwaxi 

e sin m^TTX^le Uc 


(4-57a) 


(4-57b) 


rx, O^lw.a- -jSf5T2 _ _ 

X 'e^^c ^2q sin n,irX2dx2 

Jo 

0.1wa7 , -%Lwa^, 

+ e ^ _ e e 

Jxi 

X sin niTr^ 2 ^X 2 ]djT^ 

Evaluating the integrals inside the square bracket separately: 


(4-58) 
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e e sin ni7rX2dx2 

Jo 

_ 0.1 wa- 

fX, ”07^2 

= j e cos 7i^)r2Sin n-jirx2dx2 


.Xn 1)7^2 _ _ 

-j e sin |pX2Sin n-|TTX2dx2 

J 0 ^ 


=r 

Jo 


0.1 war 
Xi Ur 2 


c ^ 

cos ^X2Sin n^TrX2dx2 
c 


O.lwa— 


^cos (n-iTT - ^)X2 - cos (n-iTr + |^)x2}dx2 


O.lwa— 

e~^^^c ,b.1wa_ . 


e I /_ _ . wa\rr /_ _ . *va\ 

_ {-g^in (n^ir + y^)x^ - (n^Tr + 


cos (n-iTT + ■Q“)xi 


+ ^^p^in (n^TT - ^)x^ - (n^TT - 5“)cos (n-jiT - ^)x-^ 


1 , wa. 

+ 2 (niTT + Uc) 
A 


1 wa 

+ 2 (n-|ir - u^) 

B 


O.lwa 

.1 . o.lwa . wa,- 

^ Uc cos (niTT - u^)x^ 
B 


+ (n^TT - )5|)sin (n^TT - g^)x^) 
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Q.lwa 1 
Ur B 


> 


0'.lwa_ 

+ iirS /0.1wa^^„ / wa.x- 

^2^ ^ Uc'^1 


+ (n^TT + J^)sin (n^TT + 5 |-)xi) 

0 . 1 wa j I 

where, ^ + (n^Tr + ]3|)^ 

n _ /0.1wa\2 . / .. wa\2 

® 

Integral Table [23] was used to evaluate the above integrals. 
Now, we want to evaluate 


wa,- 


(4-59a) 


- 0.1 wa_ wa— 

1 ~U^2 "jU^2 

e e sin niiTX 2 dx 2 

^1 


i: 


■t; 


-O.lw^ 

“Xr 


1 Ur 2 

' c wa— . j 

e cos ]]^2^''*^ niirX 2 dx 2 


-0 . 1 wa— 


rl Uc “^ 2 ] _ 

"jj_ 6 ^{cos (n-iTT - ]^)x2 - cos (n-|Tr + ^)x2^dx2 


- £ . D 
" A B 


-O.lwa- 


-e 


i)c (n^TT + g|)x^ - (n^tr + j^)cos (n^TT ^• 


-O.lwa. . 


. ”u.twa_. / Wa» — / wa\ / wa.— 

+ -jj sin (n-|TT - ^)x^ - (n-|ir - u^)cos (n^ir - 


wa> 


wa.- 
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■ 2^B^“ 


^ B T^j^os (n^TT - S|)7-, + (n^TT - jj|)sin (n^Tr - ■o^)x^> 


wa.— 


- 2 


j^-O.lwx-j 




wa^— 


wa, . 


wa,_ 


(4-59b) 


where , 

- 0.1 wa 

C = e {'O^jl^in (n^TT + u|) - (n-jiT + ^)cos (n-iir + ^)} 


- O.lwa 

D = e Uc (-0JM... 


Uc 


sin {n^TT - ^) - (niTT - ^)cos (n-iir - ^)) 


-O.lwa 


E = e (n^iT - Tg|) + (n-|TT - i 3 |)sin (n-,iT - ^)} 


-O.lwa 


F = e {-O^j^os (n^ir + ^) + (n-jU + ^)sin (n-,ir + ^)} 


Multiplying Eq. (4-59a) by and Eq. (4-59b) by and 


-0 . 1 war 


adding them, we get 


[e 


- O.lwa- O-.lw.a- , wa- 


e e ^ sin n^Trx2dx2 


O.lwa- 1 


+ eUc 


1 lJr^2 


e c s-fn n^TTX2dx2 


'1 


1 .O.lwa . /„ , wax— 1 
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. If0.1wa,.„ /- ^ wax- i 

^ ~ Uc 

-O.lwaTT O.lwa— 

1 ,1 ^x Ur ^1 ^ r Uc 1 X 

- o.lwa- Qflw a^ 

1 -» wax ~Uc j. no 1 '. 

* ■ u?'" ' 

+ (n,i, + ^)x,} 

- <"!’' - f)*l> 

- O.lwa- 0.1 wa- 

V -O.lwa Uc ""l + re"Or^l> 

+ + Fe } 


-O.lwa— 


0.1wa_ 


+ jl_{V^ ~UT^1 . Ee 
Substituting Eq. (4-59c) in Eq. (4-58) and carrying out the inte- 
gration, we get 

T r O-lwa O.lwa -. . 

"'^2AUc 2BUc 


(4-59c) 


wa 


. _ wa -O.lwa 

n, . X - g 

^ ^ yiA yiD 1-m ^ DT J 


4A 4B 

O.lwa 


^A1 B1 


+ /£_ + D_, Uc C2 Kj 

^4A 4b'® 'A1 B1 ’ 


-O.lwa 

O.lwaJ 1 x_ “DTTHl ^ FI, 

W BT^ 


LrJ J_ip 


Uc^4B 4A 
O.lwa 

+ e~W~{L- - L.}{iL - 62j 

^4A 4B^^A1 


+ j 
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+ j 


O.lwa 

2BU(.(miTr + n-iir) 

wa 




„ j. „ wa - 0-1wa 

+ J , "1^ DT ^ !2!_JJcj Uc r£l ELj 
^ 4A 4B ^B1 AV 

O.lwa 

r\ 1 1 r»n r»o 

+ J 4B^® ^A1 Bl^ 


J. ,• O.lwa rl 

J Ti t 


1 

4i 

O.lwa 


-O.lwa 


)e + — > 

4B 4A^ ^A1 Br 


+ j {fr - Ih-}6 Uc {g + 


^4A 4B^ 


'A1 B1 


(4-60) 


where 







A1 = 

(-^^-— )^ + (m-iTT - 

T 

ro 






B1 = 

(0^)2 . 

wa.2 
■ Uc^ 






Cl = 

-0.1wa„4„ f . 

— Sin (m^TT + 

wa\ 
Uc^ " 

(miTT 

+ ^)cos (niiir 


+ (m-iTT 4 


01 = 

“g^ — sin (m-iTT - 

wax 

Uc^ " 

(mi7r 

- w)cos (m-iTr ■ 

wax 

■ V 

4 (m-iTT “ 


El = 

-0.1wa^__ 

-g cos (m-iTT + 

m.) + 

lid 

(m-|7T 

+ ^)sin (m-iTT 

X W§.1 

^ Uc^ 

o.lwa 

Uc 



r"i “0*lW9 / Wd \ I f W3 \ • / Wd \ 0 • 1 W3 

= -O^osCmiTT - g|) + (m^TT - ^)sin (m^TT - jj^) + 

n - 0 • 1 wa fm ^ J. wa \ X /m X wa \ „ j ^ ^ x wa\ , 0«lw9 

G1 - -g^-^os (itiiTT + + (m-iTT + g^)sin (miir + ^) + 

C2 = 2^Msin (m^TT + ^) - (m-,Tr + ^)cos (m-jTr + ^) + (hiTt + J3|) 


02 = (m^TT - g^) - (miir - jj^)cos (m^Tr - + (ruiTr + 

no - 0*lwa„„„ /„ _ wa\ L _ W3\ . / wa\ O.lwa 

' “DT”^ ^ ^"’i^ " DH"^ ^"’l^ ■ ^1^ “ Jk^ ■ "u^ 
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Now, 


- ,0.715wb 


I 


01202 


’ I r 

• n • i 


1 fi -lyz-vil^ 


0 •'0 

Carrying out the integration we get 
1 ,0.7i5wb 


sin m 2 iTy-|Sin n 2 iryi dy-j dy 2 


I 


m2H2 - (0.715wb)2 ^ “c 


+ (nf)2Tr)(n2Tr) 


mo+no — n— 

-n + (-1) ^ 2 + e 


-0.71 5wb 


(Q. .715v<b) ^ + (^2-0)2 




(4-61) 


Then 


S"'2n,n2<^ = * Vj' 

where I^j^^^^and are given by Eq. (4-60) and Eq. (4-61) respec- 


ti vel y . 

Thus, the spectral density function of the generalized force due to 
turbulent pressure having been obtained, the generalized forces were 
simulated by using the FFT simulation technique. 

As pointed out earlier, to work out a numerical example a two 
mode approximation in the x direction and one mode approximation in 
the y direction was taken. In order to compare the results with ones 
given in ref. [l],the following numerical values taken in [1] were used 
a = 1 0 i n , , b = 20 i n . 
d = 5 in. , V = 0.3 
p = O.lh Ib/in^, E = 10^ psi 



95 


Uc = 0.65 U»» Op = 0.0056 q 
Poo = Pc = 0.00089 slugs/ft^ 
aoo = = 995 ft/sec 

ULo = 800 ft/sec, 5 = 0.157 in. 

The thickness of the plate was varied in a way to adjust the non- 
dimensional pressure, p, for the convenience of illustration. The 
damping coefficients 3-| and 82 were taken as 1 percent of the 
critical damping of the 1st and 2nd modes of linear plate, respec- 
tively. 

In order to solve the equations (4-47a) and (4-47b), 4096 
discrete values of both F](t) and F2(t) were simulated from the 
derived expression for spectral density function using the FFT simu- 
lation technique. Eqs. (4-47a) and (4-47b) were reduced into four 
first order differential equations, and a modified predictor corrector 
method was used to solve for b-|(T) and b2(T) with zero initial con- 
ditions. Fig. 16 shows the displacement time history at the center 
of the plate and the generalized forces for the first and second mode. 
The nondimensional R.M.S. response at the center of the plate versus 
the nondimensional R.M.S. pressure is plotted in Fig. 17 with and 
without the cavity effect. For the purpose of comparison, response 
corresponding to linear plate is also plotted in Fig. 17. 

Conclusion 

It is seen that the results obtained here compare very closely 
with the one given in Ref. [1] for the sub-sonic case. Results obtained 
in [1] are based on the multidimensional simulation analysis where 



Fg(t)(xio*) F|(D(xio‘) 








(T 



RMS OF NON-DIMENSIONAL PRESSURE CT- 

Fig. 17. Comparative study of plate response due to subsonic boundary layer 
pressure fluctuation by FFT method. 
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Shinozuka has taken only 100 frequency points which looks like a very 
small number considering the fact that turbulent pressure fluction 
takes place over a wide frequency range. From this point of view, 
our analysis seems more accurate because we have divided the 
spectrum into 4096 points. The total computer time to simulate 4096 
values of both F| ( t) and p 2 (t) and to obtain 4096 points for both 
b-| ( t) and b 2 (t) was only 3 minutes 30 seconds. It is needless to 
point out that the time required to simulate 4096 points for the 
generalized forces by using the trigonometric model of simulation 
will be much much more. This once again brings out the immense use- 
fulness of FFT simulation approach wherever it is applicable. 



CHAPTER V 


CONCLUSIONS AND FUTURE EXTENSIONS 

Two methods of simulation of a multivariate and mul ti correlated 
random process have been presented and these methods have been 
applied to analyze some example problems of linear and nonlinear 
random vibration. 

Comparing the two methods of simulation, we conclude that both 
the methods work very well. The only drawback is the amount of 
computer storage, which can become prohibitive if there are too many 
correlated processes to be simulated or if a very large number of 
discretized points are necessary for analysis. It is noted that FFT 
method works much faster in terms of computer time compared to the 
trigonometric model. This speed ratio is directly proportional to 
the number of series and discretized points to be simulated. For 
example, the time required to simulate the six components of wind 
turbulence, with each of the time series having 2048 discretized 
points, was twenty times less in the case of FFT method compared to 
the trigonometric approach. But FFT approach has a drawback in that 
the number of time points is equal to the number of frequency points 
and we cannot restart the simulation at any given time t. It has to 
start at t = 0. Also, the time interval is related to the frequency 
interval. The trigonometric approach does not suffer from these 
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disadvantages. Here, we can start the simulation beyond any given 
time and also the time interval is not dependent on the frequency 
discretization. Thus, we see that whereas the FFT method has a 
tremendous advantage in terms of computer time, the trigonometric 
method can be more useful in certain cases. 

As to the usefulness and importance of the simulation technique, 
we conclude that it offers probably the best means of performing the 
nonlinear response analysis of random vibration problems when the 
spectral function for the forcing function is specified. In light 
of the review study in Chapter II, we feel that the simulation approach 
and particularly the FFT method gives a much faster solution with 
lesser number of constraints. It is noted that FFT method works as 
well as the ones in Ref. [1, 6, 8] for analyzing the nonlinear re- 
sponse and this approach is much faster than that of Shinozuka 
[ 1 . 6 , 8 ]. 

The methods of simulation presented here can be extended to in- 
clude multidimensional homogeneous and nonhomogenous processes. For 
example, writing 

^ m ^ ^ ’^1 ’^2 ’ * * ’ ^ 

m Np Np 

= Z E Z E {a^ (kT,ko...kn) 

p=l ki=l k2=l kp=l 

X cos [wi|^^t + W 2 |^ 2 ^i + .... + Wpj^^Xp + o^p(kT ,k 2 ...kp)] 

X sin [wi^^t + W2,^^xi + ... + Wp^^x^ + cx^p(ki ,k 2 . . .kp) ] 

[5-1] 
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one can simulate a multivariate multidimensional process using trigo- 
nometric method. Shinuzuka [6, 8] has already extended his trigono- 
metric model to include the multidimensionality and the above proposed 
extension will not be much of a problem. But it will be really 


CAU I L I 11^ 
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the FFT approach to a multivariate-rnultidimensional 


case. It does not seem very difficult particularly when the sub- 
routine HARM can take care of up to 3-dimensional inverse Fourier 
Transfom. We feel that it should work without trouble, but it does 
need a little closer look. It might be also interesting to extend the 
FFT approach to include the nonhomogeneous process by utilizing the 
concept of evolutionary power spectrum. This will be useful in simu- 
lating trasient and siesmic response and will also lead to considerable 
time saving over Shinozuka's method [8]. 

There are lots of other areas in which the simulation approach 
can be very useful and about which we have not discussed in this 
dissertation. For example, in the area of reliability analysis and 
in the response analysis of structures when the spatial variation of 
material properties are specified. These are some of the many areas, 
the simulation technique and particularly the FFT approach can be 
extended. These extensions look very practical. 
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APPENDIX A 


1 . Random Variable and Probability Distribution Function 

Let X be any random variable which varies with statistical regu- 
larity. Since X is a random variable, any function of X is also a 
random variable. 

The first-order probability distribution function of X can be 
described by a graph as in Fig. A-1 which shows the probability den- 
sity function p(X). The hatched area gives the probability that X 
lies between a and b. It is apparent that the probability is zero 
that X actually assumes any given value x. Also we have 

j“p(X)dX=l* (A-1) 

--00 

The cumulative distribution function (CDF) defines the probability of 
occurrence of a value of X less than or equal to a given value x, i.e., 

P(x) = |\(X)dx. (A-2) 

—oo 

According to fundamental theoron of integral calculus, 

^^=p(x) (A-3) 

wherever this derivative exists. 

Properties of CDF 

(1) CDF is a non-decreasing function of x 

(2) P(x) varies between 0 and 1 as x varies from -«> to «>. thus 
p(-oo) =0 and P(<») = 1 . 

Several random variables X^ , X 2 ,...X^ are said to be jointly dis- 
tributed if they are defined as functions on the same sample description 
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space. We shall refer to p(X^, X2,...X^) as the joint probability 
density of these random variables. Thus, the probability that a^ < 
X-| < b^, a2 < X2 < b2>. . .a^ < X^ < b^ is given by 


Also, 



,...X^)dX^dX2- 


Xn)dXidX2-..dX 


..dX^. 


= 1 . 


(A-4) 


(A~5) 


The joint probability distribution of two random variables X-| and X2 
is given by 

p(XpX2) (A-6) 


and 


p(a^ < X^ < b^ 



and 82 < X2 < b2) • 
p(X^,X2)dX^dX2 


(A-7) 


Supposing that X is now a random function of spacial coordinates 
2^ or time t, the first and joint probability distribution functions 
defined above are valid for any given values of or t. That is, the 
probability that X (x or t) lies between a and b is given by 

p(X(x, or t))dX(x or t) (A-8) 

•'a 

and 


P(x) = p(X(^ or t))dX(^ or t) • 


(A-9) 


2. Mean, Mean Square, Variance and Standard Deviation 


The ensemble average 
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—00 


(A-10) 


defines the mean of X or the expected value of X. The operator E is 
used to denote this kind of average. The mean square value of X is 
given by 


E [X2] = 



(A-11) 


An important statistical parameter is the variance of X which is the 
ensemble average of the square of the deviation from the mean. Thus 
the variance of X is given by 


o2 = E [(X-yy)^] = 



(X-yx)2 p(X)dX- 


(A-12) 


An alternate expression is 

= E[X^1 - • (A-13) 

When the mean is zero, the variance is identical with mean square. 

The square root of Eq. (13), that is a, is called the standard 
deviation. The co-efficient of variation is defined by 

V = ^ (A-14) 

^X 

which tells the degree of variation. 

3. Random Process 

Let us consider random functions X(t) and Y(t) in which t 
represents any variable, say, spacial co-ordinates or time. 



no 


For a random process, we may express the correlation between 
X(t-j) and X(t 2 ) by means of the autocorrelation function 
Rx (tpt2) = E [X(t^) X(t2)l 


00 * 00 



-OO -00 


X(ti) X(t2) p(X(ti), X(t2))dX(ti) dX(t2) 


(A-15) 


The prefix auto refers to the fact that X(ti) and X(t 2 ) represents a 
product of values on the same sample at two different coordinates. 

If X(t^) and X(t 2 ) are statistically independent, we have 
p(X(t-|), X(t 2 )) = p(X(t])) p(X(t 2 )) and hence 
Rx(ti, t) = E [X(ti)] E [Xit2)l (A-16) 

Cross correlation function is defined for two samples X(t^ ) and 
Y(t 2 ) by 

Rxy (ti, t2) = E [X(ti) Y(t2)] 


00 00 


If P(X(ti), Y(t2)) dX(t]) dV(t2) 


—CO —CO 


Correlation coefficient is given by 

E [(X(t^) - E [X(ti)])_(Y(t2) - E [Y(t2)])l 
^XY = ax(ti) aY(t2) 


(A-17) 


(A-18) 


It can be shown that 

-1 < PXY £ 1 (A-19) 

and the intermediate values of pxy measure the degree of linear 
statistical dependence between X{t]) and Y(t 2 ). 



in 


4. Stationary Random Process 

A random process is said to be stationary if its probability 
distributions are invariant under a shift of the co-ordinates. In 
other words, the family of probability densities applicable at t-| also 


applies at ^ 2 ‘ particular, the first order probability density 
function p(X) becomes universal distribution independent of time. 

This implies that all the averages baseed on p(X), say, E[X] and 0 , 
are constants independent of time. 

For a stationary process, the autocorrelation or cross cor- 
relation functions become functions of the difference between t.j and 
t 2 and independent of t^ or t 2 - 

Therefore, for a stationary process , Equation (A-15) becomes 

RX(.,) = E [X(ti) X(ti+x)], T = ti - t2 (A-20) 

and Equation (A-17) is 


RxY(t) ^ [X(t-|) Y(t-| + t)], T = t] - t2* 


(A-21) 



APPENDIX B 


Fast Fourier Transform 

Restricting the limits to a finite time interval for the record 
x(t), say in the range (o, T), the finite range Fourier Transform is 
defined as 

T 

X(f, T) = J x(t)e‘^^^^^dt (B-1) 

0 

Assume now that time x(t) is sampled at N equally spaced points a 
distance h apart, then for arbitrary f, the discrete version of 
Equation (B-1) becomes 

N-1 

X(f, T) = h Z Xp Exp[-j 2tt fnh] (B-2) 

n=0 

Selecting the discrete frequency values as 

" T" k k = 0, 1, 2, ... N-1 (B-3) 


At these frequencies, the transformed values give the Fourier com- 
ponents defined by 




N-1 

Z Xp exp 
n=0 



k = 0, 1, 2 ..., N-1 


(B-4) 


where h has been included with X(fk, T) to have a scale factor of 
unity before the summation. Note that results are unique only out to 
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N 

Y. - ^ since the Nyquist cutoff frequency occurs at this point. Fast 
Fourier transform (FFT) methods are designed to compute these 
quantities, X^, and can also be used to compute the co-efficients of 
the regular Fourier series Aq and Bq in the expression 


N/2 

= x(nh) = Ao + 5 Aq cos(^^^i^) 
q=l 


N . 1 

+ 5 Bq sin(^^) . (B-5) 

q=l N 


To simplify the notation further, let 

w(u) - exp [-j . (B-6) 

Observe that w(N) = 1 and for all u and v, 

w(u + v) = w(u) w(v) • (B-7) 

Also, let 

X(k) = X|^ and x(n) = x^ (B-8) 

Then Equation (B-4) becomes 
N-1 

X(k) = E x(n) w(kn), k = 0, 1, 2, ... N-1 • (B-9) 

n=0 

Equations (B-4) and (B-9) should be studied so as to be easily 
recognized as the Fourier transform of x(n) when x(n) is expressed by 
a series of N terms. Such equations require a total of approximately 
complex multi ply-add operations to compute all of the X(k) terms 
involved. 

Fast Fourier Transform procedures are now based upon decomposing 
N into its composite (nonunity) factors, and carrying out Fourier 
transform over the smaller number of terms in each of the composite 
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factors. In particular, if N is the product of p factors such that 

N = ^r^- = r-|r2 ... rp (B-10) 

i = 1 

where the r's are all integers greater than unity, then it can be 



in an iterative fashion the sum of p terms, 

N 9 

— Fourier transform requiring 4r-| real operations each, 

N 2 

— Fourier transform requiring 4r2 real operations each, 

: (B-11) 

• 

N 2 

— Fourier transform requiring 4rp^ real operations each. 
Hence the total number of real operations becomes 


4(Nr-| + Nr 2 + Nr^ + .... + Nrp) = 4 S r 2 


i=l 

The resulting speed ratio of these FFT procedures to the standard 
method is then 


(B-12) 


speed ratio = 




_ N 


(B-13) 


4N E r,- 4 E r,- 

i=l i=l ' 


speed ratio for powers of two: 

If N = 2P, then z = 2p = 2 logo N 
i=l 

In this case, the speed ratio by equation (B-13) appears to be 

2 

speed ratio = gjjp = 
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However, a doubling of the speed can be achieved in practice by 
noting that the values for w(kn), when N is a power of 2, all turn 
out to be +1 or -1, so as to yield a higher speed ratio of the order 
of 


I a L I u 


1 _ 
4p 

_ ol 3 _ 




For example, if N = 2 '^ = 8192, Equation (B-14) gives a speed ratio 
of ^-158. 



APPENDIX C 


Expanding Eq. (3-14), with H(W|^) being a lower triangular matrix, one 
obtains 

^n "" *^iTii^n^ ^m^^n^ ‘‘‘ %n^nn’ 

^m “ l^m-il l^m2l ^ l^mml (C-1) 

From Eqs. (3-22a), (3-22b), and (3-3) and (C-1) one can see that 
m th component process is simulated from Hmp(w|() (p = 1, . . . , m) and 
the coherence between the m th component and n th component with 

n < m is generated by Hp^p(w|^) and Hnp(w|^) (p = 1, 2 n). When 

a vanishing minor is due to a zero mean square density of, say, the 
m th component process at a frequency W|^, then the elements of cross - 
spectral density associated with this component are also zero at this 
particular frequency. It is easy to show that a sufficient condition 
for this case is 

Hmn(wk) = 0, for n = 1, 2, ... m 

Hnm(wk) = 0, for n = m + l,....M 
The remaining elements of Hmn(w|^) are then determined from the sub- 
matrix obtained by deleting the m th row and m th column from the 
matrix G(w|^). 

A coherence of unity between two component processes (the m th 
and the (m + 1) th) at a frequency W|^ will also give rise to vanishing 
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principal minors. This implies a complete linear dependence at this 
frequency between the two processes and hence the existence of the 
transfer function B(w(^) such that 

G|y,+] j m+l(w|^) = |3(w) |^G|^Pfi(W|^) 


C m4-1 Atf. ^ = Q(ui\d.. (ut. \ 


(C-2) 


Therefore, at this W|^, it can be shown that G{w|^) will take the 
following form: 


^11 

®lj 

G-| , j+1 . 

••• ^iM 

^•1 


Gj » j+1 • 


sj+1 >1 • • • 

Gj+i , j 

S‘+l 

... Gj+i^M 

^M1 


G^^j+1 .. 

••• 


«ii 


^^■j 

^IM 

S'l 

S'j 

^Gjj 


^^jl 


1 G| ^Gj j .... 

^^jM 

_‘"M1 

'^Mj 





In this case, H(W|^) can be solved for a particular frequency in the 
following fashion. 

First, reduce the size of both matrices G(W|^) and H(W|^) by 
deleting the (j +1) th row and (j + 1) the column, then solve for 
the reduced H(Wj^) from the reduced G(w|(), by the method mentioned in 
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Chapter III. Finally, the deleted elements of H(W|^), as can be seen 
by direct substitution in Eq. (C-1), are given by 
Hj+1 ,p = FHjp, for p = 1 , 2, . . . j 

Hj+i,j+l=0, (C-4) 

Hp,j+1 = 0, for p = j+2, M 

If two component processes, again, say, the j th and (j + 1) th, 
are completely correlated, then 3(w) as defined for the previous case 
will reduce to a real constant 3 q independent of the frequency. This 
will lead to Dp(W|^) = 0, p = j + 1 , . . . M, for all frequencies. For 
this case, it is necessary to simulate only m - 1 component processes 
with the (j + 1 ) th removed from the system, since it can be obtained 
by 

With the proper combination of the above procedures, it is possible 
to deal with the following situations: 

(i) more than one component of G(w|^) vanishes; 

(ii) more than two components are completely coherent or 
correlated; 


(iii) the combination thereof. 
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COMPUTER PROGRAM 

This section of the report deals with the listing of some of 
the important input parameters in the computer programs included. 

On the left-hand column is the computer language (FORTRAN) symbol and 
on the right-hand side, just opposite, is the equivalent variable 
used in the derivation of the equations in Chapter III and Chapter 
IV. 


K 


Wu 

W;- 


^n 

Hji 

®mp» ^ 
Aw 


m 


^p> Op 

At 

yt) 


Listing for Program for Simulation by Trigonometric Method : 

KA, frequency index 

ALPHA(I ,KK,KA) , deterministic phase angles 

WU, upper cut-off frequency 

WL, lower cut-off frequency 

S(I,J), element of spectral density function 

HHH(I,0), elements of lower triangular matrix 

AA(I,KK,KA), BB(I,KK,KA), random Gaussian variables 

DW, discretized frequency interval 

KK,' number of time series 

VI , V2, random number generated by Gauss subroutine 
DT, time interval 

F(I,JJJ), representation of time series 
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K 

6(1, J) 
^pqk 

^pk 

^P 

h 

^pk 


Listing for Program for Simulation by FFT Method : 

KA, frequency index 

S(I,J), input spectral density function 

elements of lower triangular matrix 
Z1(K), Z2(K), independent Gaussian random numbers 
XI (K), complex random number 
DT, time interval 

DDX(I,KA), elements of complex random vector 
F(I,J), time series representation 


Of 

L 

a 

To 

P 

W] , W2 , 

A 

Sq(w) 
^mn(w) 
Tm(t) 
bn(t) 
bn(t) , 


Listing for Program for Nonlinear Response 
Analysis of String 

SIGF, standard deviation of the forcing function 

TL, length of the string 

ALPHA, a constant 

TO, initial tension in the string 

ROE, density per unit length 

W1 , W2, W3, natural frequency of vibration for 

1st, 2nd and 3rd modes respectively 

AR, area of cross section 

SO(KA), spectral density function of actual force 
S(I,J), spectral density of the generalized force 
F(I,J), generalized force representation 
AN(t), derivative of modal displacement 
BN(t), modal displacement 
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bn(t) 


^1 


32 




Uoo 

a 

b 

V 

Uc 

M 

P 

d 

u 

c 


D 

A 


ANPR(t), second derivative of the modal dis- 
pl acement 

Listing for Program for Nonlinear Panel 
Response Problem : 

BETl , percentage of critical damping for 1st mode 
of vibration 

BET2, percentage of critical damping for 2nd mode 
of vibration 

SIGP, r.m.s. pressure of boundary layer pressure 
fluctuations 

UINF, mean external flow velocity 

SA, plate length 

SB, plate width 
GNU, Poisson's ratio 
UC, convectional speed 
CM, mack number 

ROE, plate density per unit area 
SI , cavi ty depth 
SV, sound velocity 
XSI, boundary layer thickness 

ALP 

RIG, plate stiffness 
DEMC, constant 
DEM, constant 
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P SQ, dynamic pressure 

Co CONST, constant 



S iMULATIC'N OF MULT I.CORR ELATCD RANDOM PROCf'SSFS BY TRIGONMETRIC 
METHOD 

OI.MENSinN UU(6,6) ,TZ(20BG) , C16,6) ,0 ( 6 , 6 ) , AL PHA ( 6 , 6 , 6CC ) 

ltV(6,6) ,S(6,6J, U(6,6J, HH{ 6 , 6 J , HHH{ 6, 6 J , Y( 5 J , X ( 6 ) » A ( 6 ) , B { 6 ) 

1 ,AA( 6,6,600) ,BB( 6,6,600) ,W(6C0) ,WW(6) , F ( 6, 2048 ) , F I ( 20A« ) , 
iF2( 20A8) , F3(2G48) ,F4( 2043) ,FF<2G50) ,F5(2048) ,F6 (2048 ) 

COMPLEX S , HH,HHH,YYY,DD 
N=600 

WU=3. 141592 
V»L=0. 00377 
XN=N 

DW=(WU-WL)/(XN-1. ) 

D0W=S0RT(DW/6. 28318) 

0T=1. 

KA=l 

XK=KA 

DO 7 1=1,6 

DO 7KK=1,I 

ALPHA( I ,KK,KA)=0.0 

W(KA)=(XK-1.0)*0W+WL 

INPUT SPECTRAL DENSITY FUNCTION 

DO 401 1=1,6 

DO 401 J=l,6 

Cn,J) = O.G 

0(I,J)=0.0 

Z = 30.48 

CU=6. 198 

VU=0.845 

PU=(Z/18.C)**(-0.63) 

FMU = 0.03'^=Z/18.0 
US= 1.589 



H =RU=4=(US*-2)=f-'CU^'7 • 

( 1 20 . C^=F MU) * { 1 . + 1 . 5* ( W ( KA ) =1^2 / ( 60 . *2 . 0 *3 . 1 A 1 5 9<= FMU ) ) *=!■' VU ) =«''^ ( 5 . 0 

1 /(3.0^'VU)) 

ca,l)= H/G 
Z=15.2A 

BU=(Z/10.O)«#(-O.63) 

FMU=0. 03*2/18.0 
2 = 50. 

H =BU*(US**2)*CU*2 

G=( 1C1.9*FMU)* { l.+1.5*( k{ K A ) *2/ ( 5 0.9*2 .0 * 3 . 1^159*FMlJ) )**VU)**( 5.0 
1 /(3.0*VU)) 

C{4,4)= H/G 
2=30.43 
CU=8.686 
VU=0.512 

6U=(Z/18.0)**(-0.35) 

FMU=0.l*( Z/18.0)**0,58 
2=100. 

H =BU*(US**2)»CU*Z 

G=( 120.0*FMU)*( l. + 1.5*( W(KA)*2/( 60 . *2 . 0*3 . 1 41 59* FMU ) )**V'J)**(5 .0 

I /(3.o*vun 
C(2,2)= H/G 
2=15.24 

BU=U/i8.0l**(-0,35) 

FMU=0.1*t Z/18,0)**C.53 
2 = 50. 

H =BU*U'S*«? ) *CU*Z 

G=I lGi.9*FMU)* ( 1 . + 1 .5*( W(KA) *2/ ( 50.9*2.0*3. 1^1 59*FMU) )**VU)**{ 5.0 
1 /(3.C*VU)) 

C(5,5)= H/G 
2 = 1':'0. 

CCl ,3 )=(( 3. 36*1 US**2 ) *2 ) / ( 1 . 0 + 1 0 . 0* ( K A ) *2 / ( 2 . 0*3 . 114 1 5 9 * 6 0 . 0 



1 ) )«»!={5.G/3.0m/120. 
l=5C, 

C(6 ,6 ) = { (3. 36*( US *«2)*Z i / ( 1 . 0 + iC„0* ( W ( KA ) *1 / ( 2 . 0^= 3, 14 13 9^' 5 0. 9 7 

1 ))«*( 5.0/3. 0) ))/101. 94 

Z=100. 

H=(US**2)«12.54'Z 

G=-120.C*( 1 . 0*6. 0=«^W{KA)*^2/(2. 0*3. 14159^-60.0) )**(8./3. ) 

C(1 ,3)=H/G 
C ( 3 , 1 ) =C ( 1 , 3 ) 

Z=50. 

H=(US*#2)*10.0*Z 

G=-101.9*( 1 .0+6. C*W(K A) *Z/( 2.0*3. 14159*50.9) )**( 8 ./3. ) 

C(4,6)=H/G 

C(6»4)=C(4,6) 

0fLZ=50. 

DEI-F = 0.0366 
E=0.7 

C(1,4) = (C(1,1)*C(4,4) )**0.5*0.5*CEXP(-C.693*W(KA)*DELZ/i12.* 

I 3. 1416*55. 5*DELF ) J ) **C . 5* ( CGS ( W ( KA ) *0F t Z*E /55. 5 ) ) 

C(4»1)=C(1,4) 

0( 1 ,4) = (C(1. 1)*C(4,4) )**0.5*C .5*{ EXP(-0.693*W( KA) *DEL Z/'[ 2.* 

1 3. 1416*55.5*DELF ) ) ) * *C . 5* ( S I N ( W ( KA ) *DF L Z* F /5 5. 5 ) ) 

0(4,1)=-Q( 1 ,4) 

E=l,5 

PEL F=0. 0448 

C(2,5) = IC(2,2)*C(5,5) ) **0. 5*0 . 5* ( EXP ( -0 .6 93*W ( KA ) *DELZ/C 2.* 

1 3. 1416*55. 5*0tLF )) ) **0 . 5* ( COS ( W ( KA ) *DE LZ *E /5 5 . 5 ) ) 

C(5,2)=C( 2t5) 

G { 2 » 5 ) = ( C ( 2 , 2 ) *C ( 5 , 5 ) ) **0 . 5*0 . 5* ( E XP ( -0 . 6 93*W ( K A ) *OEL Z / i: 2 . * 

1 3.1416*55. 5*DELF ) ) ) * *C , 5* ( S IN( W ( KA ) *0E LZ*F /55. 5 ) ) 

0(5,2 ) =-0(2,5) 

C ( 3 , 6 ) = ( C ( 3 , 3 )*C { 6, 6 ) ) **0 . 5*0 . 5* ( EXP ( -0 . 693*W ( K A ) *OEL Z/ ( 2 . * 



I 3. 1416+35. ?=»OFLF ) I ) + +C . 5 + ( CHS ( W ( K A ) +i)F: LZ + C / 5 5 . 5 ) ) 

C{6t3)=C(3,6» 

C ( 3 , 6 ) = ( C ( 3 , 3 ) +C ( 6 , 6 ) ) + +0 . 5 +C . 5+ ( EXP ( -0 ,<S 93 + w ( K A ) +OEL Z / i[ 2 . + 

1 3. 1416 + 55.5 + DELE ) ) )*+0.5+(SIN(V'MKA) + DELZ + E/55.5) ) 

Q(o,3 )=-Q( 3,6) 

0£LF=0.04 
E = 1.0 

C(3,4) = (C(l,3)+C(4,6) ) ++0.5+0.5 + ( EXP(-0.693*W( KA ) +DELZ/( 2.+ 

1 3. 1416+55. 5+OELF ) ) ++0. 5 ) + ( COS ( W ( KA ) +DELZ + E/ 5 5 . 5 ) ) + (- 1. } 

C(4, 3)=C(3,4) 

Q(3,4)=(C(1,3)+C(4,6) )++0. 5+0.5* ( EXP ( -0 .693+W ( K A ) *DEL Z/ C 2 . * 

I 3. 1416+55.5+OELF n **.5 ) * ( S I N( W < KA ) +DELZ+E/55. 5 ) ) 

Q(4,3J=-0(3,4) 

C(1,6)=C(3,4) 

C(6, l)=C( 1,6) 

0(1,6)=Q(3,4) 

Q( 6,1) =-0(1, 6) 

DO 25 1=1,6 

DO 25 J=l,6 

C( I , J)=C< I , J)*2.0 

0U,J) = Q(I,J)*2.0 

S( I , J)=CMPLX(C( I , J) ,-0( I, J) ) 

BREAKING THE SPECTRAL MATRIX INTO THE UPPER AND LOWER TRIANGULAR 

MATRIX 

DO 29 1=1,6 

DO 29 J=l,6 

HH( I , J)=CMPLX( 0,0,0.0) 

HHH( I ,J)=CMPLX(0.0,0.0) 

HH(1,1)= CSQf-T(S( 1,1)) 

HHH(1,1)=HH(1,1) 

DO 110 J=2,6 

HH( 1, J ) = ( S( 1 , J ) /HH( 1,1)) 
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HHH( J»1 )= CCNJG(HH( 1, J) ) 

D01120 1=2,6 
CD = CMPLX( Q.0,0.0 ) 

M=l-1 

D0113 LP=1,M 
113 DD=DD+ (HH( LfM ) *C0NJG(HH( LP, I n ) 

YYY=S( I , I )-DD 
HH(I,I)= f.SQRT(YYY) 

HHH( I , I )=HH( 1,1) 

IF (I .E0.6) GO TG1121 
NN=1+1 

D0112C J=NN,6 
D0=CMPLX(0. 0,0.0) 

D0222 LP=1,KM 

222 00=00+ (HH{ LP, n*CONJG(HH(LP, J) ) ) 

HH(I,J)= ((S(I,J)-DD)/HH(I,I)) 

HHH( J, n=C0NJG(HH< I, Jl) 

1120 CONTINUE 

1121 CONTINUE 

DO 117 1=1,6 
DO 117 KK=1,6 
ALPHA! I ,KK,KA)=0.0 
AA( I ,KK,KA)=0.0 
117 BB( I ,KK,KA)=0. 0 

C RANDOM NUMBER NUMER GENERATOR SUBROUTINE GAUSS 

DO 123 KK=1,6 
CALL GAUSS! 17, 1.0, 0.0, VI) 

CALL GAUSS(29,1.C,0.0,V2) 

00 123 1=1,6 

AA(I,KK,KA)= (CABSIHHHI I,KK)) )*Vl>frODW 
BP>1 I ,KK,KA)= (CAf’:S(HHH( i,KK) ) )+-V2*DDW 



IF( A05( REAL (HHH{ I ,KK) )) .LT .0. COCCI) GO TP 2001 
ALPHA ( I ,KK , KAJ = ATAN( A [ MAG(HHH( I ,KM ) /REAL (HHH( I ,KK) ) ) 

GO TC 1213 

20C1 ALPHA ( I ,KK,KA)=1 .570B 
1233 1F(KA.GE. 20) GO TO 123 
WRITE (6, 2000) VI, V2 
123 CONTINUE 
3003 FORMAT! 21 5) 

2000 FQRMAT(2E16.8) 

IF( KA.E0.6C0) GO TO 40 
KA=KA+1 
GO TO 14 
40 1 = 1 

X X= 0 . 0 
YY=6,0 
100 TT=0.0 

DO 45 JJJ=1,2048 
45 F(I,JJJ)=0.0 

DO 46 JJJ=1,2048 
DO 47 KA=1,600 
DO 47 KK=1, I 

47 F( I ,JJJ) = F( I,JJJ)+AA(I,KK,KA)*C0S(W(KA)4TT + ALPHA( I,KK,KA) ) + 
1 BH( I,KK,KA)»SIN( W<KA)*TT+ALPHA{ I ,KK,KA) ) 

TT=TT+DT 

46 FF( JJJ)=F( I , JJJ ) 

00225 J=l,2048 
225 TZ(J)=J-1 

C PLOTTING OF THE TURBULENCE TIME HISTORIES 

FF(204<?)=16.0 
FF( 2G50)=-16.0 

CALL SCALE ( T ? , 20 46 , 2 0. 0 , TZ MI N, TX , I ) 

CALL SCALE(FF/2C5U,4.0,FFMIN,FX,1) 



CALL AXIS ( XX,YY,2H.TZ,2,20.0,C.O, TZMINtTX) 
CALL AXIS ( XX,YY,2HrF,2f4.Gf 90.0fFFMlN,FX) 
CALL PLOT (XX,YY,-3) 

CALL LINE (TZ,EF,2QA8,U 

XX=22.0 

YY=O.G 

IFir.EO.o) GO TO 111 
1 = 1+1 
GO TO 100 
111 CONTINUE 

DO 79 I=i»2C4E 
Fl( I )=F( It I J 
F2( I) = F(2, I) 

F3(I)=F(3,I) 

F4( I )=F( A, I ) 

F5( I) = F(5, I) 

79 F6(n = F(6,I) 

ViRITE(lO) FI 
WRITE(IO) F2 
WRITE(IC) F3 
WRITEIIO) FA 
WRITE! 10) F5 
WRITE! 10) F6 
CALL PLOT !O.C,C.Ct-A) 

STOP 

FNO 



SIMULATICN OF MUL T I CGRP FL ATf 0 RANDOM PPOCtSSES BY FFT METHOD 
DIMENSION UU { 6 , 6 ) , T Z ( 20 50 ) , C ( 6 , 6 ) , 0 ( 6 » 6 ) , A ( 2 048 , i , i ) 
1,V(5,6) ,SI6,6), U(6,6), HH ( 6 , 6 ) , HHh { 6 , 6 ) , FF(2050) 

It W{ 1025) , WW (6) ,P(6,204S) , Fl(2048 ) t 

1F2{ 2048) tF3I2048) , F4 ( 2043 ), F 5 ( 2048 )» F6 ( 2048 ) 

1 ,MV( 3) » INV( 1024) tZK 6) t Z2 ( 6 ) tDX ( 6 ) t ODX ( 6 , 204.8 ) , X 1( 6 ) tP( 1024) 
COMPLEX S ,HHtHHH,YYY,DD,XItGXtOOXtA 
XX=0.0 
YY=5.G 
N=1025 

WXX = 50RT(2.*204S. ) 

KU=3, 141592 
XN=N 

DW=(WU-WL)/(XN-1. ) 

DOW=SORT( DW+0.5) 

CT=1. 

KA=1 

XK=KA 

W(KA) = ( XK-1.0)’!'DW+WL 

INPUT SPECTRAL DENSITY FUNCTION 

DO 401 1= 1, 6 

DO 401 J=li6 

C( I »J)=0.0 

0(l,J)=0.0 

Z=30.4e 

CU=6.1S8 

VU=0. 945 

BU=( Z/18.0 ) *^{- 0.63 ) 

FMU = 0.03=«‘Z/1H.0 
US=1 .589 
Z=10C. 

H =RIH=( US*=<^2) 



G=< 120.0«FMU)=!= ( 1. +1 ,5*( W{ KA ) =.*^Z/( 60. *2. 0-3. 14159>:=FMU ) J =!=’!' V U )**( 5 . 0 
I /{3.0=tVU)J 
C(l,l)= H/G 
Z=15. 24 

BU=(Z/ie.O J**(-0.63) 

FMU=0.03«Z/18.0 
Z = 50. 

H =BU*{ US**2 )*CU*Z 

G = ( 1Q1.9=«=FMUJ* (l.+1.5*( W(KA)«Z/( 50.9*2. 043. 141 59--!=F HU) )**VU)**( 5.0 

1 /(3.o*vun 

C(4,4)= H/G 
Z=30.46 
CU=8.686 
VU=0.512 

6U={Z/18.0)**(-0.35) 

FMU=0.1*(Z/18.0)**0,58 

Z=100. 

H =BU*(US**2)*CU*Z 

G=( 120.0*FMU)*( l.+l .5*( W(KA)*Z/( 60, *2. 0*3. 14159*FHU) )**VU)**( 5.0 
1 /(3.C4VU)) 

C(2,2)= H/G 
Z=15.24 

BU={Z/18.0)**(-C.35) 

FMU=0.1*( Z/ 18.0) **0.58 
Z=50. 

H =BU*(US**2)*CU*Z 

G={ 101.9*rMll)*( l. + l ,5*( U(kA}*Z/< 50.9*2.0*3. 1 41 59*F MU) )**VU)** (5.0 
i /(3.0*VU)) 

C(5,5)= H/G 
Z=100. 

C(3t 3)=( ( 3.36*( US**2 )*Z ) / ( 1 . C + 1 0 . C* ( W ( K A ) *Z / (2.0*3.14159*60.0 

I ) )**(5. 0/3.0) ) )/120. 



Z = 5C. 

C(6,6)={ (3.36>!=(US«*2)<=Z ) / ( 1 . C + 1 0 . ( W ( K A ) *Z / ( 2 . 0>l'3 . 1 4159-s=5 0 . 9 7 

1 ) )'-i=*(5. 0/3.0) ) )/101. 94 

Z = 100. 

H = ( US**2) s^l2.5^'Z 

G=-120.0=>={ I .0+6.C«W(KA)=«^Z/( 2.0=4=3.14159^60.0) » =4==;= ( 5 ./3 . ) 

C( 1,3)=H/G 
C ( 3 , 1 ) =C ( 1 , 3 ) 

Z = 50. 

H=(US=4'*2) =4--10.0^--Z 

G=-101.9«( 1.0+6. C=4'W(K A) =4‘Z/( 2. 0=!'3. 14159=4=50.9) )*=!=( S./3. ) 

C(4,6)=H/G 

C(6,4)=C(4,6) 

0ELZ=50. 

0ELF=0.0366 

E=0.7 

C(1,4) = (C(1,1) =4=C(4,4) )=4'=«'0.5=s=0.5«(EXP{-0.693=«=W(KA)=4‘DELZ/(2.* 

1 3. 1416=4'55.5=('0ELF ) ) ) * =4=0 . 5=«= (COS ( W ( K A ) =«'DELZ =«=E/ 55 . 5 ) ) 

C(4,l )=C( 1,4) 

Q(l,4) = (C(lfl)=«'C(4,4) )*=«0.5=4=0.5=s=(EXP(-0.693=4=w(KA)*DELZ/(2.* 

I 3. 1416=4=55. 5=!'DELF ) ) ) . 54= ( 5 i N ( W ( KA ) *DE LZ=«=E/ 5 5 . 5 ) ) 

0(4,1 )=-G( 1,4) 

E = 1.5 

0ELF=0.044P 

C(2 ,5 )^(C( 2,2) =4=C( 5, 5 ) )’f==«=0.5’:0.5=4=( EXP ( -0 . 693=<= W ( K A ) ^OHLZ/ ( 2 .=4= 
i 3. 1416=4=55. 5=i=DELF ) ) ) =4==4=G . 5=4= ( COS ( W ( KA ) 4=DE LZ 4=F/55 . 5 ) ) 

C( 5,2)=C( 2,5) 

C ( 2 , 5 ) = ( C ( 2 , 2 J =«=C ( 5 , 5 ) ) =t'*n , 54-0 . 5=!= ( FXP ( -C . 6 93=4=W ( KA ) =4=DE I, Z/ ( 2 . * 
i 3. 1416=4=55. 5=4=nELF ) ) ) 4=4=C , 5=5= ( S I 0 ( Vv ( K A ) =4=DELZ =>'=E7 5 5 . 5 ) ) 

C(5,2)=-0(2,5) 

C(3 ,6 ) = (C( 3,3 ) =4=C ( 6, 6) ) =4=4=0. 5’4=C.5=«'( EXP{-C.693=4=W( KA) 4UELZ/( 2 .4= 
i 3. 14164=55 . 54= DELE ))) 4==4=C . 54= ( C 0 S ( W { KA ) 4=0i- LZ 4=F /5 5 . 5 ) ) 



2b 




no 


C(6,3)=C(3,6> 

Q{ 3 »6 V= t C { 3 , 3 ) \t,b) ) >!='^rj . 3<=C.S’!={ E XP (-0 .693-W { K/'-.J ^Dcll / ( 2 

i 3, 1416*55 .5«DtLF in**C.5*(SIN(W{ KA J *OELZ *E/ 55 . 5 ) ) 

Q(6,3)=-Q( 3,6) 

CELF=0.C4 
E = 1.0 

C(3,4)=(C( 1 ,3) *C ( 4,6) )**C. 5*0.5*(EXP(-0.693*W(KA)*D&LZ/(2.* 

1 3. 1416*55. 5*DELF )) **0. 5 )* (CnS(W { KA)*DELZ*E/55. 5 ))*(-!. ) 
C(4,3)=C( 3,*) 

0l3,4) = (C(l,3)*C(4,6) )**0.5*C.5*< FXP(-0.693*NtKA)*DELZ/{ 2.* 

I 3. 1416*55. 5*DFLF ) ) **.5) * ( S I Nt W ( KA ) *DELZ*E / 5 5 . 5 ) ) 

Q{4,3)=-Q( 3,4) 

Ctl,6)=C(3,4) 

C(6,n = C(l,6) 

Q(l,6)=G(3,4) 

0(6, 1)=-Q( 1,6) 

DO 25 1=1,6 

DO 25 J=l,6 

C( I,J)=C( I,J)*2.0 

0(I,J)=g(I, J)*2.0 

sn , J)=CMPLX(C( I, JJ ,-Q( I,J) ) 

BREAKING THE SPECTRAL MATRIX INTO THE UPPER AND LOWER TRIANGULAR 
MATRIX 
DO 29 1=1,6 
DO 29 J=l,6 

HH( i,J)=CMPl.X(0. 0,0.0) 

HHH( I , J)=CMPLX( O.C, 0.0) 

HH(1,1)= CSQRT(S(1,D) 

HHH( 1, 1 )=HH( 1, I ) 

CO 110 J=2,6 

HH( 1,J) = ( S( 1 ,J)/HH{ 1,1) ) 

HHH(J,1)= CCNJG(HH( 1, J) ) 


CO 

CO 



DU1120 1=2,6 
. CD=CMPLX(0.C,0.0) 

M=l-1 

D0113 LP=1,M 

113 DD=OD+ (HH( LP, I )=!-CONJG{HH(l.P,I ) ) ) 

YYY=S( I , I )-DD 

HH{ I , rl= CSGRT( YYYJ 

HHH( I , I )=HH( 1,, I ) 

IP n.FO.6) GO TO 1121 

NN= I + 1 

MM=I-1 

D01120 J=NiN,6 
CD=CMPLX(0. 0,0.0) 

00222 LP=1,MM 

222 DD=DD+ ( HH ( LP , I ) *C ON JG ( HH ( LP, J ) ) ) 

HH( I,J)= ({S(I,J)-OD)/HH(I ,D) 

HHH( J, I )=CONJG{HH(I, J) ) 

1120 CONTINUE 

1121 CONTINUE 

C RANDOM NUMBER NUMER GENERATOR SUBROUTINE GAUSS 

CO 1C<? K=l,6 

CALL GAUSS (33,0 ,707,0. C,V2 ) 

CALL GAUSS(27,0.7C7,0.C ,V3) 

Z2(K)=V2 
lO^i Z1(K)=V3 

DO 252 K=l,6 

252 yi(K)=CMPLX(Zl(K) ,Z2(K) ) 

DO 300 1=1,6 
DO=CMPLX( C.C,C.C ) 

CO 309 J=1 , I 

3C9 DO =DD + HHH(I , J) «X. 1 { J) 

DX( I )=0D 



30C DDX( I ♦KA)=DX( I ) 

IF<KA.EQ.1025) GO TC 40 

KA=KA+1 

GO TO 14 

40 00 600 1=1 t 6 

DO 6GC. N=l,1023 
K = 102.5+N 
J=10?5-N 

fOO DDX( I ,K)=CCMJG(DDX( I , J) ) 

DO 800 J=l»6 
CO 700 1=1,2048 
700 A(I ,1,1)=D0X( J,1 ) 

MV( 1)=11 
MV(2 )=0 
MV( 3)=0 

C SUBROUTINE TO PERFORM THE INVERSE FOURIER TRANSFORM 

CALL HAPM( A,MV, INV,P, 1, IFEER ) 

D0225 1=1,2048 

FFt I )=REAL(A( I, 1, 1) )/WXX 

F( J , I ) = FF( I ) 

225 TZ(I)=I-1 

C PLOTTING OF THE TURBULENCE TIME HISTORIES 

FF( 2049)=16.0 
FF( 205C)=-16.0 

CALL SCALE ( TZ , 2048 , 20 . 0 , TZM I N, TX , 1 ) 

CALL SCALE ( FF , 2 C 50 , 4 . C , FFM I N , FX , 1 1 

CALL AXIS (XX,YY,2HTZ,-2,20.,C.0,TZMIN,TX) 

CALL AXIS. ( XX,YY,2HFF,2,4.0,90.C,r-FMIN,FX) 

CALL PLCT (XX,YY,-3) 

CALL LINE (TZtFF, 2048,1 ) 

XX=22.0 



800 CONTI NUlf 

CALL PLOT (C.G,0.g,-4) 

STOP 

END 


00 



PROGRAM FOR RANDOM RESPONSE ANALYSIS OF .A NPMLINEAk STP ING USING 
FFT SIMULATIOiM HFTHOD 
niMFNSIDM W(i02A),S0( lC?5)tS(3,3) 
i,FF(2050) »F{3,20^9) , T ( 1 5 ) , Dx ( 6 ) , HOX ( 3 , 20 50 ) 

1 ,HH(4»4)»HHH(4,^f),AX(2048,l,i) ,IMV( 1024) ,Zl( 5) ,Z2(6) »Z3(6) 
i XI (6),NiV(3)tP{lC24),TZ(205C ), 

I AiPR.12048 )t A2PR( 2048) tA3»S( 2048) I 2048) ♦ S2 .Pk (2048) , 83PK(204« ) 

l#Al.MF(2C43) ,A2MF ( 2048) t A3MF (2048 ) ,B IMF ( 20-^8 ) ♦ B2MF (2048 I , P3MF( 2043 ) 

I I A 1 C R ( 2 0 4 8 ) t A 2 C R ( 2 0 4 S ) t A 3 C R ( 2 0 4 e ) , -3 1 C P ( 2 0 4 6 ) , » 2 C R ( 2 0 4 o I t B 3 C R ( 2 0 4 ft ) 
1,G1(2C48) ,G2(2048 ) ,G3(204ft) 

COMMON Al(2048) ,A2(2C48) , A3 (2048) t 83 (2048) ,B1 ( 2048 ) »b2 (2048) , 

1 AIP(2048) »A2P( 20 48) ,A3P( 2048) , B 1 P ( 2048 )» B? P ( 2043 ) ,B3P( 2.048) , 

1 FK2048) ,F2(204e)»F3(2048) 

COMMON 8ETAl,BETA2t8FTA3,Wl,W2,W3,ARtTL»T0,E 
COMPLEX AX,XI,OZ» DX,DOX 
XX=0.0 
YY=3,0 
KKK = 1 
SIG=G.2 
11 CONTINUE 

DO 239 1=1,3 
DO 239 KA=1,1025 
239 DDX( I,KA )=C MPL X ( 0 . C , 0 .0 ) 

S IGF= {SIG*42) 4-3. 14i59^n 6.0/1. 19 

N=2.-57 

XN=N 

WXX=SORT ( 2. *.029*2048 ) 

Wll=10.0*3.141.99 

WL=O.C 

DW=(WU-WL)/(XN-1. ) 

A=4 .0*3. 141 99 
TL=29.C 



AU'’HA=0.?/2f5.C 
70=100.0 
RGE = 1.0/2 5>.0 
E = 3.0-10.C<=*7 
AR=Tn/(0.05*E ) 

V.'1=SQRT( ( (1.0-3. LA150/TU*«2 )*T0/7nE ) . 

W2=SQKT( ( ( 2.0=i=3. 14159/TLl --2 )*TO/P.(>E ) 

V3 = S0RT( ( (3 .0=?3.1415.o/tl )**2 )^ro/Rnh J 

RFTAl=0.1*6.283ie/Wl 

BETA2 = 0.1=«=c.2R318/W2 

EETA3=0.1*6 .28318/U3 

KA=1 

14 XK=KA 

W(KA) =( XK-1.0)*0W+WL 

S0(KAJ = (A/( A«=i'2 + W(KA)**2n*SIGE/3.14l‘50 
JJJ = 3 

DO 15 1=1, JJJ 
DO 15 J=1,JJJ 

YJ= ( ALPHA*W(KA) )<-*2+( J#3.14159/TL )**2 
YI= ( ALPHA=4=W( KA ) M'*2 + ( 1415<^/TL )-*2 

XIJ=I*J*{ 3. 14159/TL)=!'’<'2 
lE(I.EQ.J) GOTO 18 

S( 1 ,J) = SC(KA)*X IJ*(1 .0+(-l.n )♦<•-( I+J )+EXP (-ALPHA^=W(KAJ ^TL ) 

1 *( (-1, )^=«n + i ) + ( -1. ) **( j+1) ) ) / ( Yj^'Yi ) 

GO ro 15 

C SPECTRAL DENSITY FUNCTION QP THE GF NE K AL I 7.f- D FORCt 

16 S U , J) = { ALPHA=4=W(K A) *TL+XlJ/y I=t=( l.+(-l. )*<^( I + J H E X P ( - ALPH A=i=W ( K A ) L 
1 )>:'( (-1. H-L ) + <-!. )=i=4(J+l ) n )/Yj<'Sr (K A) 

15 continue 

DO 29 1=1, JJJ 
no 29 J=l,JJJ 
HH(I,J)=O.C 



Zs' HHH( r,J) ' 

HH(1»1)= SQI^Tt S( 1, in 
HHHd , i )=HH( i, 1 ) 

DO lie J=2,JJJ 

HH( 1,J) = ( S( J ,J V/HH( 1,1) J 

110 HHH(J,1)= 

DOllZO 1=2, JJJ 
DD=0.0 

V!=I=1 . 

D0113 IP=1,M 

113 DD=00+ ( HH( LP, 1 ) 1 ) ) 

YYY= S(l,n-DD 
HH(I,I)= SQRT(YYY) 

HHH( I ,n = HHU, f J 

IF ( I ,EO. JJJ) GO TO llZl 

NN=I+1 

MM=I-1 

001120 J=NN,JJJ 
DD=0.0 

00222 LP=1,MM 

222 00=00+ (HH(LP, n*HH(l-f’,J) ) 

(( S ( 1 , J )-0'') /liHd , I) > 
HHH( J, n=HHU , J) 

1120 C3NTI!VUF 

1121 CONTINUE 

no lOP K=1,JJJ 

CALL GAUSS! 1 7 , 0 . 707 , 0. 0 , V.:l ) 

CALL GAUSS(23, J.7C7,0.C,V;:) 

Z2( K)=V2 
10<= Z1(K) = V3 

no 232 K=1,JJJ 

2 02 I ( K ) =C YH L X ( Z 1( K ) , Z 2 ( K ) ) 



no 30^ r=i,jjj 

rZ=CN*,PLX{C. 0,0.0) 

DO 309 J=1,I 

309 DZ =OZfHHH( i , J) =«=X n J) 
nx(i)=DZ 

500 DDX { I ,KA ) = DX( T ) 

IF(KA.EQ.257) GO TO AC 

KA=KA+1 

GO TO lA 

^0 DO 600 r=l,JJJ 
DO 600 N=l,1023 
K=1025+N 
J=1025-N 

600 DOX( I,K)=CONJG(DDX( I, J) ) 

00 800 J=1,JJJ 
DO 700 I=1,20A8 

700 AX(I,1,1) =ODX(J,I) - 

MV(1)=11 o 

MV(2)=0 
MV13)=0 

CALL HARM( AX,MV,INV,P,1 fl.FEER) 

SIMULATED GENERAL I ZED FORCE 
D0225 I=1,20AE 

FF( n=REAL(AX( 1,1,1) ) . 0*TI / ( TO*WXX ) 

F(J,I)=FF(I) 

225 TZn) = I-l 

FF ( 20A9 ) = 2,--. 0 
FF{ 2050)=-25.0 

CALL SCALP { TZ, 2 C'ic, 20 .O, rz>’IN.TX,l) 

CALL SCALE ( F F , 2050 , A . C , F FM I N , F X , 1 ) 

C AL L AXIS ( XX , X Y , 2HTZ , 20 . , 0 .0, TZ M I 0 , T X ) 

CALL AXIS ( XX, YY , 2HFF ,2 ,''■^.0,00.0 , FFMIN,FX ) 



o n 


CALL PLOT (XX,YY,-3) 

CALL LINF ( TZ,FF-, 2048,1 ) 

YX=22,0 

YY=O.C 

SCO CONTI NUT 

DO 801 1=1,2048 
Fl( I )=F( 1 , I ) 

F2(I)=F(2,I) 

SOI F3{I) = F{3,n 

5IMLATED GENF_F,Al I ZFO FORCE 

CALCULATION OF NONLINEAR P.ESPMSE BY NGDIFIFP/ PFEDiCTOR CORRECTOR 
1 METHOD 

C INITIAL CONDITIONS 
M=1 

H=0.025 
A1(1)=0.0 
Bl(l)=C.O 
A2I 1 )=0.0 
82{l)=0.0 
A3(l)=0.n 
B3(I)=0.0 
C STARTING VALUES 
CALL FCTIM) 

C PRFDIDTGR 

AiPk{ ?)= 2 .u4AlP( 1 ) )/3.:J 

B1PR(2)= -( 2.C*blPri ) ) /3.0 

A2PR(2)= A.n«H -(2.0 *a2P( 1 ))/3.0 

K2PP(2}= 4.04H *(2,'"*B2P(1 ) )/3.C 

A3PIU2)= 4.G<=H 2.0X=A3Pm ) /3.0 

C3PP,(2)= 4(2.C4B3P( 1 ) )/3.0 

C Nil MOOIFIFR 

Al( 2)=A1PP ( 2) 



81(-?) = 31P! (2) 

A?J 2)=A?.PR{2 ) 

B21 2)=62PP ( 2) 

A3( 2)=A3PR< 2 » 
t-3(2)=B3PP(2) 

M=2 

CALL FCT(M) 

C CORPECTOR 

A ICR( 2)=0. 12^-( p‘. C=!'A in )+ 3.0^H* ( A1 R ( 2 ) +2 . ( -A 1 P ( 1 ) ) ) 
61CR(2)=0,125*(9.0«Pl( 1 )+ 3-0*H=f ( BIP ( 2. ) <-2 .0 '-H IP ( 1 H) 
A2CP<2)=0.l25-f=<9.0*A2(l )+ 3 ( A2P( 2 ) +2 . 0-A2P { 1 ) ) ) 

B2CR( 2) = C. 125=^(9.0>fB2( U+ 3. ( B2 P ( 2 ) + 2 .0-B2P ( 1 ) ) ) 

A3CR{2)=G.125=i'(9.0«A3(l )+ 3 . ( A3P { 2 ) +2 .C <^A3P { 1) ) ) 

B3CR( 2) = 0. 125*( 9. 0=<--83( 1 )+ 3. ( B3P ( 2 ) +2 . C ^83 P ( 1 ) ) ) 

C FINAL VALUES 

A1(2)=A1CP( 2H-9.0 « ( A1 PR < 2 )- AlCR ( 2 ) ) /12 1 • 0 

81(2) = BlCR(2)+9.0 <={ 31 PR ( 2 )- BICR ( 2 )) /1 21 . C 

A2( 2) = A2CR(2) + 9.0 ♦( A2PR { 2 ) - A2C.R ( 2) )/l21.0 

32(2)=B2CR(2J+9,C *( B2PR( 2 I-B2CR ( 2 ) ) / 121. 2 

A3( 2I=A3CR(2J+9,G =«'( A3 PR ( 2 )- A3CP ( 2 1 1/121.0 

H3( 2) = B3CR( 2) + 9.0 «( 33PR ( 2)-B3CR ( 2) ) / I 2 1. 0 

C THIRD POINT 
N=2 

WR I Tt ( 6 , 66 ) ni ( -M ) » P2 ( N ) » P3 ( N ) , F 1 ( M ) , F2( M ) , F 3 { M J 
CALL FCT(M) 

A1PR(3)=4.0 *H^(2.0«A1P(2)-A1P(1))/3.0 

B IP P, ( 3 ) =4 . 0 i 2 . 0=^B 1 P ( 2 ) -0 IP m ) / 3 . C 

A2PR( 3) = 4. 2 s^H=i=( 2.0’f'-A2P(Z )-A2 p (I ) )/3.C 

HZPR { 3) = A.O 2.0*B2P( 2)-B2P.( 1 ) )/3.^. 

A3PR(3)=4.0 2.0*A3P(2 )-A3P (1 ) )/3.0 

fi3PR(3)=A.n 2.0<=B3P( 2)-o3P( 1 ) >/3.C 

r Mf ni F IBR 



.Uf't H ( 3 ) =A i PR ( 3 ) - 1 1 2 . 0 A 1 PR ( 2 ) - A 1 CR { 2 ) ) / 1 2 1 .0 

8lMF(3):=RlPPr3)-l 12.0 *(B1PR(2)-B1CR(2) )/121,0 

A2MF(3) = A2PR( 3 )-l 12.C * ( A2PR ( 2 ) - A2CR ( 2 ) ) / i 2 1 . 3 

(32i'^Ff 3) = 82PR(3J -1 12.0 « ( B2PR ( 2 ) - B2CR ( 2 ) »/ 1 2 1 . 0 

A3MF(3) = A3PR(3)-112.0 { A3PP ( 2 ) - A3CR I 2 ) ) / 1 2 1 .0 

B3MF( 3)=a3PP(3)-112..0 ’SM R3PR ( 2 ) -B3CR ( 2 1) / I 2 I .0 

AU3> = A1MP(3) 

Bl( 3)=aiMF( 3) 

A2( 3) = A2i'1F(3) 

62( 3)=B2MF(31 
A3( 3)=A3MF<3) 

B3( 3)=B3MF( 3) 

M=3 

CALL FCT(M) 

A1CR(3) = 0.125=P(9.0*A1{2 ) + 3.0*H*{ A1P( 3) + 2.n#AlPJ2)-AlP( 1) ) » 

81CR(3) = 0,125#(9.0’!=B1(2)+3.0’>=H«(81P(3)+2.0*B1P(2)-B1P( 1) )) ^ 

A2CR(3)=0.125’^(9.0«A2(2) + 3.0*H*{ A2P(3)+2.0*A2P(2)~A2P11J ) ) w 

B2CR (31=0, 125^(9. C*B2 (2 1+3 . 0=«=H« ( B2P (3 1+2.0*B2P( 2 )-B2P( 1)1) 

A3CR(3)=0. 125* ( 9. 0* A3 ( 2 1 + 3, C*H* ( A3P ( 3 ) +2 . 0*A3P ( 2 1 -A3P ( 1 ) ) ) 
e3CR(3) = 9ol25*(9.0*B3(2)+3.0<=H«(B3P(3)+2.0*B3P(2)-B3P{ 1)1) 

C FINAL VALUES 

AI(3)=A1CR( 31+9.0 *( A IPR C 3 1 - AlCR ( 3 ) ) / 1 2 1 . 0 

Bl( 3)=B1CP ( 3 1+9.0 *(P1PR(3)-B1CR( 3) 1/121.0 

A2(3l = a2CF(3U9.0 «(A2PR(31-AZC.R(311/121.0 

82( 3)=B2CR(31+9.0 *( 32PR ( 3 1- B2CP ( 3 1 1 / 1 2 1 . C 

A3( 3)=A3rR( 31+<^^.C * ( A 3 PR ( 3 1 - A3CR (3 1 1/121.0 

E3(3}-=R3CR131 + 9./; *(83FRl 3)-B3CR(51 )/12i.O 

C FOURTH POINT 
M = 3 

WR 1 1 E ( 6 , 66 1 N, 8 1 ( ,3 1 , 82 ( ^ , B3 ( P: 1 » F 1 { N ) » F 2 ( N 1 , F 3 ( M 1 

CALL rr.KNi 

AlOR ( R}=R.O*H «( 2.0*AiP(3)-AlP(2l-2.0*AlP( U 1/3.C- 



B 1 P R ( 4 ) = . 0 ’I - H =? { 2 . 0 * B 1 P ( 3 ) - B 1 P ( 2 > - 2 . U B 1 ! > ( 1 ) ) / 3 . G 

A 2 u H ( 4 ) = 4 . 0 V 1 1 =iM 2 . 0 * A 2 P ( 3 ) - A 2 P ( 2 J - 2 . 0 * A 2 ( 1 ) ) / 3 . 0 

B2PR ( ^ ) =A. 0-AH *(2.0* B2P ( 3 ) -P2P ( 2 ) - 2 . G*B2 P ( I ) ) / 3 . 0 

A3PR( 4J-=4.0*H *( 2.0*A3P(3)-A3P(2 )-2.0*A3P( 1) )/3.0 

B3PR(A)=A. 0*H *( 2.0*B3P(3)-H3P<2 )-2.0*i33P( 1 ) )/3.0 

66 FORMAK I5t6F16.6) 

C MOniFIER 

A I M F ( A ) = A 1 P R ( 4 ) - 1 1 2 . 0 * ( A 1 P R ( 3 ) - A I C R ( 2 ) ) / 1 2 1 . 0 

81MF( A) = 31PR( A J-1 12 .0 * I B1 PR ( 3 )- B1 CR ( 3 ) ) / 1 2 1 .C 

A2MF(A)=A2PR(A)-11Z.G *(A2PR (3)-A2CR(3) ) / 1 2 1 . 0 

B2MF(A> = 82PR(A)-112.0 * ( R2PR < 3 )-R2CR ( 3 ) ) /I 2 1 . 0 

A3MF ( A) =A 3PR ( A ) - 1 12 . 0 * ( A3P« ( 3 ) - A3CP. < 3 ) ) / 1 2 1 . 0 

R3N'F(A) = B3PP(A)-1 12.0 * ( B3PR { 3 )- B3CR ( 3 ) ) / 1 2 1 . 0 

Al( A)=A1MF ( A) 

Bl( A)=B1MF( A) 

A2( A)=A2MF(A) 

62( A)=82MF( A) 

A3( A)=A3MF(A) 

B3( A)=B3MF(A) 

M=A 

CALL FCT(M) 

A1CR{ A)=0.125*( 9oO*Al(3 )- Ai ( 1 ) +3 . 0*H* ( A IP (A) + 2.G*AiP(3)-AlP(2))) 
01CR( A)=C. 125*(9,C*BI (3 ) - Rl ( 1 ) +3 .C*H* ( f3 1 P ( A ) + 2 . 0 *B 1 P ( 3 ) -B IP ( 2 ) ) ) 
A2CR(A)=0. 125*( 9.0*A2(3)-A2( 1 ) +3.0*H* ( A2P ( A ) + 2 . 0* A2 P ( 3 ) - A2 P ( 2 )) ) 
A3CR( A) = 0. 125*(9.0*A3{3)-A3( 1 )+3.0*H*( A3!M A)+2.0*A3P( 3)-A3P( 2) ) ) 
B3CR( A)^0. 125*( . 0*B 3 ( 3 )- B3 ( 1 ) +3. 0*H* ( b3P ( A ) + 2 . 0* B3 P ( 3 I - B3 P ( 2 i) J 

R2CR ( A»=0.12 5*( 9.0*B2(3 )~B2 ( 1 ) +3 . 0*H* ( B2P ( J + 2 . 0 *B2 P ( 3 i -B 2P ( 2) ) ) 
C FINAL VALUES 

A1(4)=A1CR( A)-9.0 * ( A 1 PP ( A )- A iCP ( 4 ) ) / 1 ;> i . 0 

B 1 ( A) =n IC R ( 4) -9 . 0 *( B1 PP ( ^ )-BlC2 ( '- ) ) / 1 2 1 . 2 

A2( A) = A2CR( A)-9,(: *( A2Pf A? CR ( A ) ) / 1 2 1 . ( 

A3 ( A ) =A3CR ( A )-s . 0 *( A3PR ( ‘A )- ABTR ( ^ ) ) / 1. 2 1 . 



B3(4)“H3CR{ A)-9.0 =<=( B3PR ( )-83CK ( 4) ) / I 7. 1 .0 

F32( 4) = B2CR(4)-9.0 h2 F R ( ^ - B2CR ( ^ ) ) / 1 2 1 . 0 

DO 139 I = , 2045 
M=I 

CALL FCT( M) 

A1PR( M+i)=A 1( N'-3) +4.0«H =.M 2. C*A1 P ( R ) - A 1 P ( M- 1 ) + 2 . A 1 P ( M-2 ) )/3.0 

BIPR ( H+1) = B1 < M- 3 ) +4 .0«H ^ ( 2 . O^Bl P { M ) -3 iP ( H- L ) + 2 . 04'r. 1 p ( 2 ) ) /3,C 

A2PR (M+ 1)=A.2( F;-3 )+4. 0-H t-' ( 2. C ’i' A2 P ( M ) - A 2P ( M-1. ) + 2 . 0 * A2 P ( M-2 ) ) / 3.0 

B2PR( M+1 ) = B2( ^-3 ) +4.0#H ( 2 . 0*02 P ( M ) -B2P ( M- I ) + 2 . 0«!;;2 P ( ,R-2) )/3.0 

A3PR(M+l)=A3{N-3)+4.0-^H * ( 2 . 04A? P ( M ) - A3 P (M-1 ) + 2. A3 P ( M-? ) )/3,0 

B3PR{M+l) = B3{M-3)+4.0=f=H »( 2. C*B3P ( H ) -33P ( M- 1 ) + 2 . 0*33 P ( M-2 > )/3.0 

C MODIFIER 

A1MF{MH) = A1PR(M+1)-112.0 * ( A1 PR ( M ) - A ICR ( M ) ) /121.0 

B1MF{M+1) = B1PR( M+U-112.0 * ( B1 PR ( M ) -B1 CR ( M ) ) / 12 1 . 0 

A2MF(M+l)=A2PR(M+i)-H2.0 # ( A2PR ( M) - A2CR ( M ) ) / 1 2 1 . 0 

B2MF(M-H) = B2PR(M + 1)-112.0 *{ B2PR( M )-B2C.P ( M ) ) /121 .0 ^ 

A3MF(M+1)=A3PR(M+1)-112.C « ( A3PR ( M ) - A?. CR ( M ) ) / 12 1 . C 

R3MF(M+l) = B3PR(M+l)-il2.0 « ( R3PR ( M )-B3CR ( M ) ) / 121 . 0 

AH M+1) = A1MF(M+1) 

8l( M-H) = B1MF(M+1 ) 

A2( M+1»=A2MF(m+i) 

R2( M+1)=B2MF(M+1 ) 

A3( M+l)=A3MF{M+l ) 

B3< M+1 )=33MF( M+1 ) 

■■ ■D=M + 1 

CALL FCTIN) 

A1CR(M+1) = 0, 125* ( 9. 0*A1 {MJ-A 1 (M-2 ) +3.U*H* ( A IP ( M+1 )+2 , 0*A1P (M)-AIP 
i (M-im 

D1CR(M+1)=0. 125* ( 9.0*B1 ( M ) -R1 ( R-2 > + 3 . C *H* ( 0 I P ( M+ 1 ) +2 .0*3 IP ( M ) - 5 IP 
1 (M-l))) 

A2CR ( M+1) ^0, 12 3* { 9. 0*A2 (M)-A2 (M-2 ) + 3.M:*H* ( A2P ( M+l ) + 2. (M ) - A2P • 

1 (M-1))) 



B2CR(M + 1 )=0.:12'5- ( 9.0«B2 ( '■1 )- B2 ( M- 2 ) +3 . ( I? ?P ( 4- 1 ) + 2 . j=i-B 2 R ( B ) - R 2 P 

1 (M-im 

A3CR ( M+ 1 ) = 0 . 1 2 5=!= ( ? . 0 *A3 ( M ) - A 3 ( H- 2 ) + 3 . 0=f=h* ( A 3P ( 1 ) + 2 . 0* A 3 P ( K ) - A3 P 

1 ( M-in ) 

B 3CR ( M+ 1 ) = Q.. 1 2 5- ( 9.0*83 ( M ) - B3 ( M- 2 ) +3 . 0*H* ( S 3 P ( ’'■1+ 1 ) +2 . 0* B3 P ( M ) - B3 P 

1 

C FINAL VALUES 

Al( M+1 )=A1C.R( M + 1 J-9.0 *( AlPR (.3 + 1 )-AlCF (N+j. ) ) /121 .0 

ei( 3+l)=81CP ( w + l )-9.0 *(81PR(M + 1)-[UCR (M + i) )/l2i.O 

A2( M+1) = A2CP.( N+1 )-9.0 * ( A2PR ( M+ I )- A2CP ( M+ 1 ) ) / 1 2 1 . 0 

B2( M+L)=B2CR( R+1 )-9.0 *( B2PR (M+1 )-82CR ( m + i ) ) /121 .0 

A3( M+1) = A3CP (M+U -9.0 * ( A3PR ( M+1 ) - A30P. ( M+1 ) ) /12L .0 

B3(M+i)=83f.R(N+l)-9.0 #( B3PR (M+i )-83CR ( 3 + U ) / 121 .0 

139 CONTINUE 

00 llOl I=l,2CA6 
1101 AK n=Rl( I )-B3( I ) 

00 201 1=1,2046 

201 TZ(I)=I 

C PLOTTING THE RESPONSE FOR DIFFERENT MODES 

31(2047) = !. 0 
Bi( 2048)=-1.0 

CALL SCALE ( TZ , 20^6 , 20. 0 , TZ M I N, TX , 1 ) 

CALL SCALE ( B 1 , 20 48 ,4 . 0 , 3 1 M 1 N , 6 X , I ) 

CALL AXIS( XX,YY,2HTZ,-2 ,20 . 0 , 0 .0 , TZ M I N , TX ) 

CALL AXIS (XX,YYf2HRl,2,4.0,90.0,BlMIN,3X) 

CALL PL0T(XX,YY,-3) 

CALL LINE(T?,B1,20^6,1) 

XX=22.0 

YY=0.0 

DO 202 1=1,2046 

202 TZ(I)=I 

B2( 2047)=1.0 



B?.{ 20AR)-^-l. 0 

CALL SCALB ( r2,2'.:iA6t20.0,TZ'N!lN,TX,l) 

CALL SCALf (i32,2OA3,4.C,02MIM,BX,l) 

CALL AX I S ( X X » Y Y ♦ 2 HTZ f -2 , 2 0. 0 , 0.0 , TZ K-I N , IX ) 

CALL AXIS ( XXtYY,2H22.?,4.0,90.0,82MIfNi,PX) 

CALL PLOT ( XX, YY,-3) 

CALL LINEC TZ,C2,2046, 1) 

XX=22.0 

YY=O.C 

DO 203 1=1,2046 
203 TZ(I)=I 

B3( 2047)=1.C 
B3( 2048)=-! .0 

CALL SCALE ( TZ , 20 46 ,20. 0, TZM I N, TX , 1 ) 

CALL SCALE ( 3 3 , 2048 , 4. 0 , S3M J N t BX , 1 ) 

CALL AXIS<XX,yY,2HTZ,-2,20.0,0.0,TZMIN,TX) 

CALL AXIS ( XX,YY,2HB3,2,4.0,90.0,B3MIN,BX) 

CALL PLGT(XX,YY,-3J 
CALL LINE (TZ,B3, 2046,1) 

00 205 [=1,2046 
20 5 TZU) = I 

C PLOTTIMG OF TE!£ RESPONSE AT THE CEMTEP OF THE STRING OUE TO ALL 

1 THE 3 NOnES GOMBINEO TOGETHER 

A1(2047)=1.C 
Ai( 2043)=-! .0 

CALL SCALE [ T Z , 20 46 , 20. 0 , TZM I N, TX , I ) 

CALL SCALE ( Ai , 20 4P , 4 ,G , AIXII N, AX , 1 ) 

CALL AxrS( XX,YV, 2HTZ,-2,2C.G,0.G,TZNIN, TX ) 

CALL AXIS ( XX, YV,2HAl, 2,4. 0,90.0, AIHIN, AX) 

CALL PLOT ( XX,YY,-3) 

CALL L INF (TZ,Al ,2046, 1) 

CALL PLOT (0.0, 0.0 , -4) 



i 


STOP 

f?NO 


SUSfiOUTINE: USFO TO CALCULATE THE FIF-ST SCCONF) DERIVAIIVES 

OF THE RFSPONSF DOF TO DIFFFRFNT OODES 
SUBROUTINE FCT(M) 

COMMON AK2048 ) , 42(2048 )t A3(2C43 ) ,t3l 204 B) »F'1( 204S ) ♦E-2( 204.B) , 

1 A1P( 20AR) , A2P( 2C4S) , A3P( 204» ) ,B1 P( 204 8 ) , '32 P( 204.8 ) , R.3P( 2043 ) , 

1 FK2048) ,t^2(2048)»F3(204B) 

COMMON BFTAl»BFT.A2,BFTA3,wl,W2»iv3,AR,TLtTCs E 

A IP { Ml=-2.0*0. 628 31 8’^ A I ( M ) - ( I . 0+ ( ( AP v h ) / (4. j*T;; J ) * ( 3 . 14 1 5F/T L) *=."2 
1 *( B 1 ( M N) =«’-*2+9.G’>B3 (M )*^=2 ) )'-^{ 6 . 2 S3 1 8 4^2 ) 'M-. U .m )+F 1 ( M) 

61P(M)=A1(M) 

A2P (M)=-2 .0*0.628318*42 (M)-( i.0 + ( ( A1'4E )/( 4 .04T0n 4( 2. 1415 9/TL ) *=1=2 
1 *(B1(M)**2+4.0*B2( M)**2+9.0*B3 (K J442) ) *( 12 .566 36**2) *B2 < M ) ♦ F 2 ( M ) 
62P(M)=A2(MJ 

A3P (M)=-2. 0*0.6283 18* A3 ( M )-( i.0+( ( AR*f; ) / ( A . 0 *T 0 ) ) *( 3 . 141 54'/TL ) **2 
1 *(B1(M)**2+4.G*B2(M)**2+9.0*B3(M)**2))*( 18,84954**2) *H3 ( M ) +F3( M ) 
B3P (M)=A3(M) 

RETURN 

END 



C PROGRAM FCR NONLINEAR PANEL RESPONSE DUE TO TURBULENI PRESSURE 

1 FLUCTUATIONS FOR SUBSONIC BOUNDARY LAYER 

CIHENSIGN C(2t2),Q(2,2) ,W(2048I, F ( 2t 3000 ) , HH{ 2 , 2 ) tHHH( 2 , 2 I » 

1 AX(4100, 1,1) , INV(2048) »Z1(3) ,Z2(3) tS(2»2) »TZ ( 1505) , 

1 XI(3), MV(3),nDX(2,41G0),DX(3),FF(300C), 

1 A1PR(3000),A2PR(3000),B1PP(30001,B2PR(3000), 

1 AIMFOOOG) ,A2MF(3000),B1MF(3000),B2MF(3COO), 

1 AlCROOOO) ,A2CR(3000),B1CR(3000) ,B2CR(300C) 

COMMON AK3000) ,A2( 3000), BK 3000), 82 (3000), 

1 A1P(3000),A2P(3000),B1P(3000),B2P(3000), 

1 FU3000),F2(3000),8FT1 ,8ET2, 

1 C10,C11,C12,C21,C22,C33,C44,C20,C66,C55 
COMPLEXS, DYY, HH,HHH,YYV,DD,DDX, AX , XI , DX 

C COEFFICIENTS IN THE GOVERNING EQUATIONS 

BETl=0.5 
BET2-1.7 
KOUNT=l 

3532 SIGP=7500«0 
K0UT=1 

399 DO 401 1=1,2 
DO 401 J=l,2 
C(I,J)=0,0 

401 Q(I,J)=O.C 

DO 239 1=1,2 
DO 239 KA=l,204o 
239 CDX( I,KA)=CMPLX(0.0,0.0) 

WU= 1000. 0^3. 1^159 
WU= 2000.0*3.14159 
WL= 0.0 

H4=56. 0*128. 0*8.9 *12. 0*0.91/ ( ( 10 .0**7 ) *l44 , 0*S 1 GP ) 

H=H4**0.25 

TH1=3.14159*C.5 



TH2=3. 14159*0.5 

UINF=800.0 

XSI=0. 157/12. C 

SA=in.o 

SB=20.0 

Sl=5.0 

GNU=0,3 

RO=0. l*H/386.4 
N=2049 
XX=0.0 
YY=3.0 

UC=0.65*UINF*12.0 
CM=800. 0/995. C 
ALP=0.5 

PI4=(3.14159)**4 

ALPS=ALP**2 

GNUS=GNU**2 

ROC=O.OOC09 

SV=995.0 

XN=N 

DW={MU-kL)/(XN-1.0) 

KA=1 

E =10. 0**7 

RIG=E*(H**3)/(12.0*( l-GNU^^Z) ) 

SQ=8. 9*32.0 

0EMC= R0C*SV**2*SA*SB**2/(RIG*144.0) 

DEM=2.0*SQ*SA*SB**2/( RIG*144.0) 

C10=( (ALP+l.O/ALP )**2)*PI4 
C20=( ( ALP+4.0/ALP )**2)*PI4 

C11=0.75*PIA*( ( 1.0-GNUS)*( ALPS+1.0/ ALPS) +2.0* (2.0 *GNU+ALPS+1.0/ 
1 ALPS)) 

C12=0.75*P1 4*11 1 .0-GNUS )*(4.0*t ALPS +1 .07 ALPS ) 1 . 0*AL PS 



1 /( 1.0+A.0*ALPS)<'*?+ALPS/(9.0+4.0*ALPS)=«'*2)+2.0«=( 5.0^GNU + ALPS + 

1 4.U/ALPS J J 
C21=C12 

C22=0.75=«'PI4*( ( 1.0-GNUS )«(ALPS*16. 0/ALPS) +2. 0*( 8. 0*GNU+AL.PS+ 16.0 
1 /ALPS)) 

C33=DEMC#(-SA/(PI4*S1) ) 

C44=-DEM*0.666/CM 
C66=DEM*0. 6666/CM 

C55 = S0RT(RIG) /(4.«CM*S6*UINF<'SORT(RO)<'12.0) ♦DEM 
SD2= SQRT(RIG/R0)/{ SA*SB) 

C0NS1= ( SIGP**2)«XSI«SD2/< {3.14159**2)*3«136»UINF ) 
DELT=SD2/(1000.0) 

DELT=SD2/(2000.0) 

WXX=S0RT( DELT*4096.0) 

WRITE(6.466) RO»SO, DEMC , DEM ♦ SD2 , RI G ,0ELT,WXX 

466 E0RMATaX,8F13.6) 

WRITE (6, 196) C10,Cll»C12,C21,C20,C22,C33fC44,C55,C66 
196 fORMAK 10F12.5) 

SIG=H 

INPUT SPECTRAL DENSITY FUNCTIONS OF THE GENERALIZED FORCES 
14 XK=KA 

W(KA) = (XK-1.0H'DW+WL 
C0NS2=3.7*EXP(-2.0*W(KA) ♦XSI/UINF) 

C0NS3=0.8*6XP( -.47»W(KA)*XSI/UINF ) 

C0NS4=3.4*EXP(-8. 0*W(KA) ♦XSI/UINF) 

CON ST^CONSK't CONS 2+C0NS3-C0NS4) 

Ml=l 
M2 = l 
Nl=l 
N2=l 

CQ=0.i«W(KA) =>10.0/00 
ES=EXP(-0.1*W(KA) *10.0/00 



EX=EXP(0.i*W(KA)*10.0/UC) 

AG1,'N1*3.1A159+W(KA)«10.0/UC 

AG2=N1«3. 14159-W(KA)*10.0/UC 

AC=CQ**2+AG1*«2 

BC=CQ«*2+AG2=«‘*2 

SAG1=SIN( AGl) 

CAGl=COS( AGl) 

SAG2=SIN( AG2) 

CAG2=CQS( AG2» 

C CONSTANTS ON PAGE 91 ' 

CC=ES*(-CQ*SAG1-AG1«CAG1) 
DC=ES*(-CQ*SAG2-AG2*CAG2) 
EC=eS*(~CG*CAG2+AG2«SAG2) 
FC=ES*(-CQ<'CAG1+AG1*SAGU 
C CONSTANTS ON PAGE '«3 

AA1=M1*3. 14159+W(KA)=C'10.0/UC 
88l = Ml*3.1A159-W(KA)<‘10.0/UC 
SA=SIN(AA1) 

SB=SIN(BBH 

CA=COSC AAIJ 

CB=COS(BBU 

CC1=-CQ«SA-AA1«CA+AAI 

0D1=-CQ=»'SB-BB1 AC6+BB1 

CC2=CQ*SA-AA1*CA+AA1 

DD2=CQ*SB-DBI*C8+BB1 

EE1=-CQ*CA+AA1«SA+C0 

FFl=-CQ«CB+BB14'SB+CO 

GG1= CQ«CA+AAl«SA-CO 

GG2= CQ*CB+Bei«SB-CO 

AC1=C0*#2+AA1««2 

BC1 = C0*<'2 + RB1««2 

ACC l=( 1A.30«W(KA) /UC )=«‘*2+(N2*3.1A159)*=<=2 



BCC1=( 14.30*W{ KA) /UC1*«2+{M2«3.14159)=S=*2 

C( 1 1 U = 0.5*CQ/AC + 0. 5^C0/SC + (0.25*AG1/AC+0.2 5*AG2/BC )*ES« 

1 (CCl/ACl+DOl/BC l) + {0.25=^CC/AC+0.2 5*DC/BC)^EX«{i:C2/ACl + DD2/8Cl ) 

1 +CQ*(0.25/BC-0.25/AC)*ES=«‘(EEl/ACl + FFl/BCl) 

1 +(0.25«FC/AC-0.25*EC/BC)*EX*(GGl/ACl-GG2/eCl) 

1 +( 14.30*Vi(KA) /UC+( (3.1A159)*«2*<H2*N2)/BCC1)*{ 1 . 0 + ( - 1 . 0 ) ( M2+N2 ) 

1 +EXP(-14.30*W(KA )/UC) ♦( ( -1 . 0) <=* ( M2+ U + ( -1 . 0 ) ♦« ( N2+1 ) ) ) )/ACCl 
C(I,1)=C{ 1,U*C0NST 

Q( l»l)= (C.25*AG1/AC+0.25*AG2/BCI*ES«‘(EE1/AC1 + FF1/BC1 ) 

1 +(0.25«CC/AC+0.25*0C/BCI«EX*(GG1/AC1+GG2/BC1 ) 

1 +CQ«(0.25/BC-0.25/AC)*ES*(CCl/ACl+DDl/BCl) 

1 +(0.25#FC/AC-0.25*EC/BC)#EX«(CC2/AC1+DD2/BC1) 

0(lf 1)=0(1,1 )^CCNST 

Ml=l 

M2=l 

Nl=2 

N2=l 

AG1=NI*3.14159+W(KA)*10.0/UC 
AG2=N1*3.14159-W< KA)*10.0/UC 
AC=CQ<'*2+AG1**2 
BC=CQ=«'*2+AG2**2 
SAG1=SIN(AG1) 

CAG1=C0S(AG1) 

SAG2-SIN( AG2) 

CAG2=C0S( AG2) 

C CONSTANTS ON PAGE 

CC=ES=S'(-CC*SAG1-AG1=!'CAG1) 

DC=ES*(-CQ<'SAG2-AG2*CAG2) 

EC=ES=«'(-CG<'CAG2 + AG2«SAG2) 

FC=ES-«'(-CQ*CAGl + AGl*SAGl ) 

C(l,2)= ( 0.25^AGl/AC+0.25«AG2/BC)«ES* 

1 ( CCl/ACl + DDl/BCl ) + <0.25*CC/AC+0.25<'DC/BC)=!‘EX#(CC2/ACl+D02/BCl ) 



1 +CQ*(0.25/BC-0.25/AC)*ES<'(FFl/ACl + FFl/6Cl) 

1 +(0.25*FC/AC-0.25«FC/BC)^EX«(GG1/ACL-GG2/BCI ) 

1 +1 14.30*W( KA)/lJC + ( (3.14159)*«2*(M2*N2)/BCC1)*( 1 . 0 + { - 1 . 0 ) ** ( M2+N2 ) 
1 +EXP(-1A.30*W( KA) /UC) ’!=( ( - 1 . 0) *♦ ( M2 + 1 ) + ( -1 . 0 ) ** ( N2 + 1 ) ) ))/ACCl 
C(1,2)=C( 1,2)=«'C0NST 

Q( 1 ,2)=CQ/( AC<'(M1«3.14159+N1*3.14159) )-CQ/ ( BC*( Ml *3 . 1 4159>N 1* 

1 3. 14199) ) 

1 +(0.25=»=AG1/AC+0.25«AG2/BC)*ES«( EE1/AC1+FF1/8C1) 

1 +( 0.25*CC/AC+0.2 5«0C/BC)=J'EX*(GG1/AC1+GG2/BC1 ) 

1 +CQ*(0.25/BC-0.25/AC )=«=ES4=(C.Cl/ACl + Dni/BCl) 

1 +(0.25*FC/AC-0.25*EC/BC)*EX«(CC2/AC1+DD2/BC1) 

Q(1,2)=Q(1,2)«C0NST 
C(2,1)=C( 1,2) 

Q(2,l)=-0(1,2) 

Ml=2 

M2=l 

Nl=2 

N2=l 

AA1=M1«3.14159+W(KA)«10.0/UC 
B81=M1*3.14159-W( KA)«10.0/UC 
SA=SIN( AAl) 

SB=SIN(BB1) 

CA=C0S(AA1) 

CB=C0S(BB1) 

CC1=-CQ*SA-AA1*CA+AA1 
D01=-CQ«SB-B81 *CB+BBl 
CC2=CQ*SA-AAI<=CA+AA1 
DD2=CQ«SB-8B1*CB+BB1 
EE1=-CQ*CA4-AA1«SA+C0 
FF1=-C0«C8+361<'S8+C0 
GG1= C0*CA+AA1*SA-CQ 
GG2= CQX'CB+PBl’f'Se-CQ 



f 


AC1 = C0'^«2+AA1*«2 
BCI=CQ<=*2 + RB1»«2 

C(2t2)=0.5'«'CQ/AC + 0.5=i=CQ/8C+ (C.25«'AG1/AC+0.25*AG2/BC)*ES* 

1 (CC1/AC1 + D01/BC1) + (0.25<'CC/AC + 0.25*DC/BC)*EX«(CC2/ACI+DD2/BC1 ) 

1 +CQ*(C.25/BC-0.25/AC)«ES«(EEl/ACl+FFl/BCl) 

1 +{0.25«FC/AC-0.25«EC/RC)«EX«(GG1/AC1-GG2/BC1 ) 

1 +( 14.30^W( KA) /UC+l (3,14159)«*2*(M2*N2)/BCCl)«(1.0+(-1.0J*«(M2+N2) 
1 +eXP(-14.30*W(KA)/UC) ♦( (-1.0)«*(M2+l) + (-1.0)*=«'(N2+l ) ) ) )/ACCl 
C(2»2)=C( 2,2)*CONST 

Q(2,2)= (0.25*AG1/AC + 0.25*AG2/BC)*ES’^(£E1/AC1 + FF1/BCI ) 

1 -K0.25^CC/AC+0.25*DC/BC)*£X«(GG1/AC1+GG2/BC1 ) 

1 ♦CQ<'(0.25/BC-0.25/AC)*ES*(CCl/ACl+D0l/8Cl) 

I +(0,25*FC/AC-0.25*FC/BC)*EX#(CC2/AC1+D02/BC1) 

Q(2,2)=Q(2,2)*C0MST 
DO 25 1=1,2 
DO 25 J=l,2 

25 S(I,J)=CMPLX(C(I,J),-Q( I,J)» 

C CALCULATIGNS OF THE GENERALIZED FORCES STARTS 

JJJ = 2 

DO 29 1=1, JJJ 
DO 29 J=1,JJJ 
HH( I,J)=CPPLX(0.0,0.0) 

29 HHHl I , J)=CKPLX{0.G,C,0) 

HH{1,1)= CSQPT(S(1,D) 

HHH(1,1)=HH(1,1) 

DO 110 J=2,JJJ 

= (S( l,in 

110 HHH(J,1)= CGNJG(HH( 1, J) ) 

D01120 1=2, JJJ 
DD=CKPLX(0. 0,0.0) 

M=I-1 

D0113 LP=1,M 



113 DD=DD+ (HH( LP, I ) ^CONJG(HH(LP, 1 ) ) ) 

YYY= S( It I )-0D 
HH( I , I )= CSCRT(YYY) 

HHH( I I, I) 

IF ( I .FP.JJJ) GO TO 1121 

NN=1+1 

MM=I-1 

D01120 J=NN,JJJ 
D0=CMPLX(0. 0,0.0) 

00222 LP=1,MM 

222 00=00+ (HH(LP,I)<'CONJG<HH(LP,J))) 

HH(i,j)= nsn ,j)-DD)/HHn , I) ) 

HHH( J , I ) =CON JG ( HH ( I , J ) ) 

1120 COMTINUE 

1121 CONTINUE 

DO 109 K=1,JJJ 

CALL GAUSSn7,0.707,0.C,V3) 

CALL GAUSS(23,0. 707,0. 0,V2) 

Z2(K)=V2 

109 Z1(K)=V3 

00 252 K=1,JJJ 

252 XUK)=CNPLX(Z1 (K) ,Z2(K) ) 

DO 300 1=1, JJJ 
DD=CMPLX(0. 0,0.0) 

DO 309 J=1,I 

309 CD =DO+HHH( 1,J)*XI(J) 

DX(I)=OD 

300 D0X( I ,KA) = OX( I) 

IF( KA.EQ.2049) GO TO 40 

KA=KA+1 

GO TO 14 



on 600 N= 1,2047 

K=2049+N 

j=204Q-N 

6C0 DDX( l,K) = CaNJG(OOX( I, jn 
DO 800 J=1,JJJ 
DO 700 1=1,4096 
700 AX( 1,1,1) =DDX(J,I) 

MV( 1 j=l2 
MV( 21=0 
MVt3)=0 

CALL HARM( AX,MV,INV,P,1,IFEER) 

D0225 1=1,3000 

FF( I)=REAL{AX(I,l,l))/WXX 

F(J,I)=FF<n 
225 CONTINUE 

800 CONTINUE 

DO 801 1=1,3000 
Fl(I)=F(l,I) 

801 F2(I)=F(2,II 

C SOLUTION OF THE GOVERNING DIFFERENTIAL EQUATIONS BY MODIFIED 

1 PREDICTOR CORRECTOR METHOD STARTS 

C INITIAL CONDITIONS 
M=1 

H=OELT 
All 11=0.0 

Bim=c.o 
A2I1)=0.0 
B2(ll=0.0 
C STARTING VALUES 
CALL FCT(M) 

C PREDIDTOR 

A1PR(2)= 4.0«H <^(2 .C*AIR< 1 ) )/3.3 



r 


B1PR(2)= ^.C^H «(2.C*B1P( 1) )/3.0 

A2PR(2)= 4,C«H ?.0*A2P( 1) )/3.0 

B2PR(2J= 4.0*H ♦(2.0*32P( l))/3.0 

C NO MODIFIER - 

A1(2)=AIPR(2) 

B1(2)=81PR<2) 

A2( 2)=A2PR(2) 

B2( 2)=82PR( 2) 

M=2 

CALL FCT(M) 

C CORRECTOR 

A1CR(2)=0.125«(9.0*A1(1)+ 3.0*H*( AlP ( 2 )+2.0«AlP( 1)) ) 
B1CR(2)=0.125»(9.0«B1( 1 )+ 3.0«H* ( B1P(2)+2,0*81P( 1 ) ) ) 
A2CR(2) = 0.12 5^{9.0«A2( H+ 3. 0*H»( A2P( 2 >+2 .0*A2P ( IM ) 
82CR(2) = 0.125*(g,0«B2( U+ 3.0*H* ( B2P ( 2) +2.0=«'82P ( I )) J 
C FINAL VALUES 

Al( 2)=A1CR( 2K9.0 ♦( A1PR(2)-A1CR (2) ) /121.0 

Bl( 2 )=B1CR( 21+9.0 ♦( 81PR ( 2 )-BlCR ( 2 ) ) /121.0 

A2( 2)=A2CR(2)+9.0 *( A2PR ( 2 1 -A2CR ( 2 ) ) /1 2 1.0 

B2(2)=B2CR{2)+9.0 «( B2PR ( 2 )-B2CR ( 2 ) ) /12 l.O 

C THIRD POINT 
M=2 

WRITE (6f66) M, Bl( M),B2tM), F1(M),F2(M) 

CALL FCT(MI 

AlPRt3)=A,C *H«(2.0«A1P(2)-A1P( 1 l)/3.0 

B1PR(3» = A.C <=H«(2.0+B1P(21-BIP( 1) )/3,0 

A2PR(3) = 4.0 «H*(2.n=«'A2P{2)-A2P(l) )/3.0 

B2PR<3)=4.0 2.0*82P(2)-B2P(1) )/3.0 

C MnOIFIER 

A1MF(3)=A1PR(3)-112.0 ♦ ( AlPR ( 2 )-AlCR ( 2 ) )/121.0 

B1MF(3)=B1PR(3)-I12.n ♦(B1PR(2)-B1CR(2))/121.0 

A2MF(3) = A2PR(3)-112.0 ♦ ( A2PR ( 2 I- A2CR ( 2 > ) / 1 2 1 .0 



f 


B2MF(3i = B2PR(3)-112.n { B2PP. ( 2 ) -B2CR ( 2 ) ) /1 2 1 . 0 

AU3)=AlMFt3) 

BI(3)=B1MF<3) 

A2(3)=A2MF(3) 

82(3)=B2WF(3) 

M=3 

CALL FCT(M) 

A1CR( 3)=0.125«(9.0*A1(2 ) + 3.C*H*< A1P( 3 )+2.0*AlP(2 )-AlP( 1) il ) 
B1CR( 3) = 0. 125*( 9.G*Bl(2)-»-3.0«H«(BlP(3)+2.0*BlP{2)-BlP( 1 ) ) ) 
A2CR( 3)=0. 125«(9.0*A2 (2 H-3.0«H*( A2P(3)+2.0*A2P(2)-A2P( !») ) 
B2CR( 31 = 0. 125*( 9,0>>B2(2) + 3 .C«H=f(B2P(3) + 2.0*B2P(2)-32P( 1) ) ) 
C FINAL VALUES 

All 3)=AlCR(3)+9.0 *( A1 PR ( 3 )- A1 CR ( 3 ) ) /1 2 1 . 0 

Bl( 3)=81CR(3)+9.0 *( B1PR( 3 )-BlCR( 3) )/121.0 

A2(3)=A2CR(3)*9.0 ♦( A2PR (3 ) - A2CR ( 3) )/ 121.0 

B2(3)=02CR(3)+9.O ♦( B2PRI 3 )-B2CR ( 3 )) /121. 0 

C FOURTH POINT 
M=3 

CALL FCTIM) 

A1PR(4)=A.0*H ♦(2.0*A1P(3)-A1P(2)-2.0*AIP( 1) )/3.0 

B1PR(4)=A.C«H «(2.0«81P(3)-B1P(2)-2.0*B1P( l))/3.0 

A2PR( A)=A.O*H *( 2.0*A2PI3)-A2P(2)-2.0*A2PI 1) )/3.0 

B2PR( A)=A.O=FH <=( 2.C*B2P(3)-B2P(2 )-2.0*B2P( 1) 1/3.0 

66 FORMAT! I5fAF16. 6) 

C MODIFIER 

A1MF(A) = A1PRIAJ-1 12.0 « ( A1 PR I 3 ) -A ICR ( 3 ) ) /12 1 . 0 

BImFI A)=B1PRI<^)-112.0 *(81PR(3)-B1CRI3))/121.0 

A2MF( A)=A2PR(4)-112.0 ( A2PR ( 3 ) -A2CR ( 3 ) ) / 1 21 . 0 

B2MF( 4) = 82PR(4)-112.0 * ( 82PR ( 3 )- B2CR ( 3 ) ) /121 .0 

Al( 4)=A1MF( 4) 

R1(4)=B1MF(4) 

A2( 4)=A2MF( A) 



62{ ''+) = B2MF( A) 

M=4 

CALL FCT(M) 

A1CR( 4)=0. 125«( g.G^'Al (3 )-Al( I ) +3 . ( A 1 P ( 4 ) +2 . 0 =«=Ai P ( 3 ) -A IP ( 2 ) ) ) 
R1CR(4) = 0. 125*(9.0=frRl(3)-31( 1)+3.0=«'H<'(BIP(4)+2.0*B1P(3)-BIP(2) ) ) 
A2CR( 4) = 0. 125=t (9.C*A2 ( 3 )-A2 ( I )+3 .0«H* < A2 P ( 4 ) +2. 0«A2P ( 3 )- A2P { 2 ) ) ) 
B2CR( 4)=0. 125*( 9. 0*B2{3)-B2( 1 )+3. ( B2P ( 4 ) +2 . 0* B2 P ( 3 ) -82P ( 2 ) ) ) 


C FINAL VALUES 

Al( 4)=A1CR( A)-9.0 
Bit 4)=BlCR(4)-9.0 
A2(4)=A2CR(4)-9.0 
B2(4)=82CR(4)-9.0 


<--(AlPR(4)-AlCR(4) )/121.0 
*(31PR(4)-B1CR(A) )/12l.O 
*( A2PR(4)-A2CR(4) )/121.0 
♦<82PRI4I-B2CR(4) ) /121.0 


DO 139 1=4,2999 

IF( I.GT,2500) WRITE(6,56) FI ( M) ,F2 ( M ), Bl ( M) 
56 F0RMAT(3F16.6) 


M=I 

CALL FCTIMJ 

AiPR(M-H)=Al(M-3)+4.0*H 
BlPR(M+l) = Bl(M-3)+4.0<'H 
A2PR(K+l)=A2(M-3)+4.0*H 
B2PR(M+1)=B2(M-3»+4.0«H 
C MODIFIER 

A1MF(M+1)=A1PR(M+1)-112.0 
B IMF (M+1)=B1PR(M+1) -112.0 
A2MF(M+1)=A2PR(M+1»- 112.0 
B2MF(M+1)=Q2PP (H+1) -112.0 
All M+1)=A1MF(M+1 ) 

Bl( M+1)=81MF(M+1 ) 

A2( M+1)=A2MF(M+1J 
B2(M+1)=B2MF(M+1) 

N=M + 1 


*(2.0*A1P(M)-A1P(M-1)+2.0«A1P( M-2) )/3.0 
♦(2.04BlP(M)-BlP(M-l)+2.0*BlP(M-2) )/3.0 
«(2.0*A2P(M)-A2P( M-1)+2.0*A2P(M-2) ) /3.0 
♦ (2.0*B2P(M)-B2P(M-n + 2.0#B2P(M-2) 1/3.0 

AlPR(M)-AlCR(Mn/121.0 
=«'(BlPR(M)-BlCP(Mn/121.0 
A2PR(M)-A2CR(wn /121.0 
B2PR(M»-82CR(M}) /121.0 


CALL FCT(N) 



» 


T 


A1CR( M+i)=0.12 5«( P.G*A1 ( M ) - A 1 { M-2 ) +3 { A IP ( .M + 1. ) +2 . O’^'A 1 P J M ) - A IP 
I 

B1CR(M+1)=0.125=C'( g.O«Bl (M)-Bl(M-2J+3.0=i'H*(BlP( M+1 )+2.0*BlP(M)-BlP 
1 (H-1))) 

A2CR(M+1) = 0.12 5*( 0*A2 ( M)- A2 ( M-2 J + B.O’i'H* ( A2 R ( M+ 1 ) + 2 . D<= A2P (M ) - A2P 

1 (M-im 

82CR( M+l )=0.125* ( 9.0*B2 ( M )-B2 (M-2 ) +3 . { B2P ( M+ 1 ) +2 . 0«B2P ( M ) - B2P 

1 (M-im 

Al(M+l)=AlCR(H+l)-9.0 « { A1 PR (M+ I ) -A ICR ( M+1 ) )/121 .0 

Bl{M+l) = BlCR(M+l)-9.0 * ( BIPR ( M+ 1 1-BlCR ( M+1 ) > /121.0 

A2< M+1)=A2CR ( M+U-9.0 « ( A2PR ( M+ 1 ) -A2CR ( M + 1 ) )/ 121.0 

B2{M+l)=B2CR(M+l)-9.0 « ( B2PR (M+1 )-B2CR ( M+1 ) ) /I 21.0 

IF(BKM) .GT.IOOO. ) GO TO 732 
139 CONTINUE 
332 11=2047 

USUM=0.0 

DO 1 1=1500,3000 
1 USUM=USUM+B1( 

XII=3000.C-1500.0 
RMS = SQRT(USUM/Xin 
WRITE(6,96) SIG,RMS 
96 F0RMAT(2F16.8) 

555 FORMAT! 16) 

732 WRITE<6,555) M 

C PLOTTING OF THE SIMULATED GENERALIZED FORCES AND THE RESPONSE 

1 AT THE CENTER OF THE PLATE 

XX=0.0 
YY=3,C 

DO 2271 1=1,500 

11=1+2499 

TZ(n=I 

Fl( I ) = F1( T I ) 



F2{ I ) = F2( I I ) 

2271 Bl(I)=8ini) 

81(501) = 2.0 
BK 502) = -2.0 

CALL SCALE ( TZ , 500 ,2 0. 0 , TZ M IN, TX , 1 ) 

CALL SCALE { B 1 , 502 , 4. 0 , B IM I N, BX , 1 ) 

CALL AXIS(XX,YY,2HTZ,-2,20.0,0.0,TZMIN,TX) 
CALL AXIS (XX,YY,2HB1,2,4.0,90.0,B1M1N,BX) 
CALL PL0T(XX,YYt-3) 

CALL LINe(TZ,81,5CC,l) 

XX=22.0 

YY=0.0 

DO 3271 1=1,500 
3271 TZ(I)=I 

F1I501)=800.0 
Fl( 502) =-800.0 

CALL SCALE ( TZ , 500, 20.00, TZK IN, TX , 1 ) 

CALL SCALE ( FI , 502,4 .00 , FIM IN,FX, 1 ) 

CALL AXIS(XX,YY,2HTZ,-2 ,20.G ,0 .0 , TZ MIN, TX ) 
CALL AXIS ( XX,YY, 2HF1,2 ,4.0,90.0,F1MIN,FX) 
CALL PLCT(XX,YY,-3I 
CALL LINEITZ, FI ,500,1) 

XX=22.0 

YY=0.0 

DO 4271 1=1,500 
4271 TZ(I)=I 

F2( 501)=800.0 
F2( 502)=-800.0 

CALL SCALE ( T Z , 500 , 2 0 . 00, TZM I N, TX , 1 ) 

CALL SCALE I F 2 , 50 2, 4. OC , F2M1 N ,F X , 1 ) 

CALL AXIS! XX,YY,2HT Z,-2,20. 0,0.0, TZM I N,TX) 
CALL AXIS (XX,YY,2HF2,2,4.0,90.0,F2MIN,FX) 



CALL PLOT (XX, YY, -3) 
CALL LlKE(TZ,L2,500il > 
CALL PL0T(0,0,0.0,-A) 
STOP 
END 



SUBROUTINE USFO TO FIND THE FIRST AND SFiCONn DERIVATIVES 
FOR SOLVING THE GOVERNING DIFF. EQUATIONS 
SUBROUTINE FCT(m) 

COMMON Ai(50CC)tA2(300CI,BlI30CO),82(3000)t 
1 A1PI3G00) ,A2P(3G00I ,B1P{30CO),B2P(3000) , 

1 Fl(3000),F2(300G).BETl,8ET2» 

1 C10,Cll,C12»C21,C22»C33tC4A,C20,C66,C55 
A1P(M)= -RET 1«A1(M)-(C10+C1 1*81 ( M)«=S^2+C1 2*B2( M) =>=!=2 ) *B1 (M)+4,.0«{ 
1 F1(M)*-C33*4.0*B1(M) ) 

81P(M)=A1(M) 

A2P(M)=-BET2*A2(M)-(C20+C21*B1IM)**2+C22*B2(M)**2)*R2(H)+4.0*( 

1 F2(M)) 

B2P(M)=A2(MJ 

RETURN 

END 
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